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ABSTRACT 


Monte  Carlo  procedures  designated  as  LITE-I  and  LITE-II  were 
developed  to  study  the  transport  of  light  through  the  earth's  atmos¬ 
phere  under  various  environmental  conditions.  LITE-I  treats  monochroma¬ 
tic  light  emitted  from  a  point  source ,  and  LITE-II  treats  monochroma¬ 
tic  plane  sources  of  light.  The  codes  have  been  written  In  both  ALGOL 
for  the  Burroughs  B-5000  and  FORTRAN  II  for  other  computers.  The  codes 
are  sufficiently  flexible  to  treat  multiple  scattering  In  an  atmosphere 
In  which  air  density  and  aerosol  size  distribution  vary  Independently 
and  arbitrarily  with  altitude.  Provision  for  treating  ground  and  cloud 
reflection  with  an  albedo  method  is  also  available  In  the  codes. 

The  codes  have  been  verified  through  comparisons  with  other  cal¬ 
culations  of  light  transport  In  the  atmosphere.  Utilization  Instruc¬ 
tions,  Input  data  formats,  sample  problems,  and  the  ALGOL  listings  of 
the  codes  are  given  to  aid  those  who  wish  to  utilize  the  codes. 


ill 


FOREWORD 


The  authors  wish  to  express  their  appreciation  to  Henrietta 
Hendrickson  and  Hemma  Francis  of  Oak  Ridge  National  Lsbcratory  who 
aided  in  the  checkout  and  running  of  test  problems  on  the  FORTRAN 
versions  of  the  LITE  codes.  They  alao  wish  to  acknowledge  the 
assistance  of  Leon  Leskowitz,  of  the  U.  S.  irmy  Electronics  Laboratory, 
Fort  Monmouth,  New  Jersey,  for  his  assistance  in  translating  the 
FORTRAN  codes  to  the  ALGOL  language  and  his  many  helpful  suggestions 
during  the  checkout  of  the  ALGOL  versions  of  the  LITE  codes.  The 
work  described  in  this  report  was  carried  out  under  the  technical 
monitorship  of  Dr.  R.  W.  Fenn  of  the  Atmospheric  Science  Laboratory, 
USAECOM,  Fort  Monmouth,  New  Jersey. 


iv 


TABLE  OF  CONTENTS 


Abstract  Hi 

Foreword  iv 

List  of  Figures  vii 

List  of  Tables  ix 

I.  Introduction  1 

II.  Preliminary  Investigations  6 

III •  Methods  19 

3.1  Description  of  LITE-I  Code  22 

3.1.1  Selection  of  Source  Angles  22 

3.1.2  Optical  Path  Lengths  25 

3.1.3  Atmosphere  Description  28 

3.1.4  Geometry  29 

3. 1.4.1  SEARCH  Routine  31 

3. 1.4. 2  Distance  To  a  Boundary  31 

3. 1.4.3  Region  Importance  Numbers  33 

3.1.5  Statistical  Estimation  34 

3. 1.5.1  Detector  Angle  of  Scatter  36 

3. 1.5. 2  Probability  of  Scattering  Toward  a 

Detector  37 

3. 1.5. 3  Contribution  From  a  Single  Collision  37 

3.1.6  Light  Scattering  38 

3. 1.6.1  Rayleigh  Scattering  39 

3. 1.6. 2  Mle  Scattering  41 

3. 1.6. 3  Direction  After  Collision  41 


v 


3.1.7  Ground  and  Cloud  Reflections  44 

3. 1.7.1  Reflection  Into  a  Detector  45 

3.1. 7.2  Reflection  Angle  47 

3.1.8  Printout  of  Scattered  Intensities  47 

3.1.9  Direct  Intensity  Calculation  49 

3.2  Description  of  LITE-II  Code  50 

IV.  Comparison  of  Lite  Code  Results  With  Other  Data  54 

4.1  LITE-I  Code  Results  54 

4.2  LITE-II  Code  Results  69 

References  86 


vi 


LIST  OF  FIGURES 


Figure  Page 

1.  Dependence  of  the  Scattering-to- total  Cross  Section  Ratio  on 

the  Albedo  8 

2.  Dependence  of  the  Scattered  Flux  at  4-4  Mean-Free-Paths  from 

a  Point  Isotropic  Light  Source  on  the  Maximum  Number  of 
Collisions  Allowed  Per  History  9 

3.  Comparison  of  Calculations  of  Photon  Flux  Times  Distance 

Squared  for  Infinite  Medium  and  Air-Ground  Geometries  13 

4.  Variation  of  the  Total  Flux  After  a  Given  Number  of  Colli¬ 
sions  with  Penetration  Distance  15 

5.  Light  Intensity ®Received  by  a  Flat  Plate  Receiver  as  a  Func¬ 
tion  of  Slant  Range:  Ground  Albedo  -  0,5  16 

6.  Comparisons  of  Calculations  of  Light  Transport  m  which  the 

Ground  Albedo  was  Taken  to  be  0.5  and  0.0  18 

7.  Comparison  of  Cumulative  Distribution  of  Angles  Sampled  Using 

Rejection  Technique  With  Cumulative  Rayleigh  Phase  Function  42 

8.  Comparison  of  LITE-I  Code  and  K74  Results  for  an  Infinite 

Medium  55 

9.  Comparison  of  Hand-Calculated  and  LITE-1  Code  Single  Scattered 

Intensities:  t  *  0.5,  8  *  0°.  Albedo  -  0,0  58 

c  * 

10.  Scattered  Intensity  Transmitted  Through  a  0.13  Mean-Free-Path 
Thick  Variable  Density  Rayleigh  Atmosphere:  Normal  Incident 

Line  Beam  Source  60 

11.  Scattered  Intensity  Transmitted  Through  a  0.5  Mean-Free-Path 
Thick  Rayleigh  Atmosphere:  Normal  Incident  Line  Beam  Source , 

Albedo  *  0,8  62 

12.  Scattered  Intensity  Transmitted  Through  a  0.5  Mean-Free-Path 
Thick  Rayleigh  Atmosphere:  Normal  Incident  Line.  Beam  Source, 

Albedo  ■  0.0  63 

2 

13.  4irR  times  Intensity  as  a  Function  of  Distance  from  a  Point 

Isotropic  Source  in  an  Infinite  Medium:  Scattenng-to-Total 
Cross  Section  Ratio  =  0.5  65 

14.  Intensity  Versus  Distance  from  a  Plane  Isotropic  Source  in 
an  Infinite  Medium:  Scatterlng-co-Total  Cross  Section  Ratio 

-  0.5  67 

vi  i 


Figure  Pap,e 

15.  Scattered  Intensity  as  a  Function  of  Distance  from  Plane 

Parallel  Monodirectional  Source  in  an  Infinite  Medium:  6  * 

0%  Scattering-to-Total  Cross  Section  Ratio  ■  0.5  0  68 

16.  Angular  Distribution  of  the  Scattered  Radiation  Transmitted 
Through  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  “ 

0°,  Ground  Albedo  -  0.0  °  72 

17.  Angular  Distribution  of  the  Scattered  Radiation  Transmitted 

Through  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  * 

0*,  Ground  Albedo  -  0.25  0  73 

18.  Angular  Distribution  of  the  Scattered  Radiation  Transmitted 

Through  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  = 

0°,  Ground  Albedo  *  0.8  °  74 

19.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5 

Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  -  0°,  Ground 

Albedo  -  0.0  °  75 

20.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5 

Mean-Free-Path  Thick  Rayleigh  Atmosphere;  0  ■  0°,  Ground 

Albedo  -  0.25  °  76 

21.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5 

Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  *  0°,  Ground 

Albedo  -0.8  0  77 

22.  Angular  Distribution  of  the  Scattered  Radiation  Transmitted 

Through  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  » 

cos"^0.6.  Ground  Albedo  =0.0  °  78 

23.  Angular  Distribution  of  the  Scattered  Radiation  Transmitted 

Through  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere,  0  = 

cos“^0.6.  Ground  Albedo  =0.25  °  79 

24.  Angular  Distribution  of  the  Scattered  Radiation  Transmitted 

Through  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere,  0  ■ 

coa“l0.6.  Ground  Albedo  «0.8  °  80 

25.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5 

Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  =  cos“^0.6. 

Ground  Albedo  “0.0  °  81 

26.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5 

Mean-Free-Path  Thick  Rayleigh  Atmosphere:  3  «  cos“^0.6. 

Ground  Albedo  “0.25  °  82 

27.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5 

Mean-Free-Path  Thick  Rayleigh  Atmosphere:  0  “  cos“^-0.6, 

Ground  Albedo  ■  C.8  °  83 

28.  Spatial  Distribution  of  the  Radiation  Intensity  as  a  Func¬ 

tion  Depth  into  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmo¬ 
sphere:  Normal  Incident  Plane  Parallel  Source  85 


vii:. 


LIST  OF  TABLES 


Table  Page 

I.  Dependence  of  Various  Parameters  on  Maximum  Number  of 

Collisions  Allowed  Per  History  in  L05  Problems  11 

II.  Scattered  Light  Transmitted  Through  Planetary  Atmospheres 

for  a  Normal  Incident  Plane  Farallel  Broad  Beam  Source  59 

III.  Scattered  Radiation  Emerging  from  Upper  and  Lower  Surfaces 
of  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere  Due  to  a 
Normal  Incident  Plane  Source  70 

IV.  Scattered  Radiation  Emerging  from  Upper  and  Lower  Surfaces  of 

a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere  Due  to  Plane 

Source  Incident  at  0  »  cos"l0.6  70 

o 


1 


I.  INTRODUCTION 

In  order  to  predict  the  degree  of  thermal  radiation  damage  that 
may  occur  from  nuclear  detonations,  the  effects  of  various  atmospheric 
and  terrain  conditions  upon  the  propagation  of  light  energy  must  be 
evaluated.  During  recent  years,  considerable  emphasis  has  been  placed 
on  both  theoretical  and  experimental  aspects  of  light  transport  in 
the  atmosphere.  Typical  of  the  experimental  data  are  observations  re¬ 
ported  by  Canter  and  Petriw  (Reference  1) .  They  measured  the  scat¬ 
tered  light  intensities  from  a  point  source  above  partially  and 
totally  snow-covered  terrain  in  which  the  surface  albedo  varied  from 
approximately  0.2  tc  about  1.0.  Some  of  thexr  observations  were 
made  during  cloud  cover,  providing  data  on  the  "duct  effect"  or 
energy  trapping  situations.  Curcic,  Kenes trick  and  Cosden  (Refer¬ 
ences  2,  3,  4  and  5)  made  a  thorough  study  of  light  scattered  by  the 
atmospheric  aerosol.  Many  of  their  observations  were  conducted  over 
the  water  surface  using  a  collimated  source  and  receiver  arrangement 
that  reduced  the  detection  of  any  surface  reflected  light. 

The  U„  S.  Naval  Radiological  Defense  Laboratory  has  also  con¬ 
ducted  experimental  studies  (References  6,  7,  8  and  9)  on  the  trans¬ 
mission  of  light  emitted  by  a  4i?  source  through  the  atmosphere.  The 
results  of  their  measurements  give  information  on  the  scattered-to- 
direct  intensity  as  a  function  of  both  source-receiver  range  and 
receiver  angle  for  the  source  and  receiver  located  near  the  air- 
ground  interface.  The  USNRDL  measurements  give  data  on  the  scatter¬ 
ing  of  light  in  the  atmosphere  for  atmospheric  conditions  ranging 
from  very  clear  (visibility  15+  miles)  and  for  various  conditions  of 


cloud  cover. 
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In  general,  these  experimental  data  indicate  that  the  transport 
of  light  within  the  earth's  atmosphere  is  highly  dependent  upon  both 
the  meteorological  and  terrain  conditions.  This  dependence  is  so 
great  that  it  is  difficult  to  apply  measurements  made  at  one  time 
and  place  to  predict  the  scattered  light  intensity  at  some  other 
time  or  place.  Experimental  investigations  of  the  intensity  of 
light  at  a  given  receiver  point  as  a  function  of  a  single  parameter, 
svch  as  the  surface  albedo,  are  hampered  by  the  inability  of  the  in¬ 
vestigator  to  control  all  of  the  other  environmental  parameters. 

Several  of  the  experimental  data  for  point  sources  in  the  atmos¬ 
phere  have  been  used  to  derive  semi-empirical  relations  which  will 
allow  prediction  of  the  light  intensity  as  a  function  of  distance 
from  an  uncollimated  source.  Typical  of  these  semi-empirical  rela¬ 
tions  are  those  derived  by  Stewart-Curio  (Reference  12) ,  and  Gibbons 
(Reference  8).  All  of  these  semi-empirical  relations  correlate  the 
uncollimated  transmission  with  the  extinction  coefficient  and  the 
path  length.  To  account  for  ground  albedo  or  cloud-cover  reflections, 
some  correction  factor  is  usually  applied  to  these  relations  <■  At 
the  present  time,  there  is  no  general  agreement  on  the  magnitude  of 
the  correction  due  to  ground  and  cloud  albedo  and  the  dependence  of 
the  correction  factor  on  the  source-receiver  distance,  or  the  wave 
length  of  the  source  radiation t  It  should  also  be  noted  that  there 
is  considerable  disagreement  between  the  results  obtained  from  these 
semi-empirical  relations  at  distances  larger  than  a  few  kilometers. 

Much  of  the  theoretical  work  published  on  light  scattering  has 
obvious  limitation  for,  in  many  cases,  the  calculational  methods 
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developed  are  beset  with  difficulties  when  extensions  are  made  to 
other  than  simple  models.  The  Boltzman  transpot  equation,  which 
rigorously  defines  the  transport  of  light,  is  generally  written  in 
intractable  form,  and  many  fidelity -reducing  procedures  have  been 
applied  to  the  transport  equation  in  order  to  obtain  solutions. 

These  include  replacing  an  inhomogeneous  medium  with  a  homogeneous 
medium,  the  use  of  Isotropic  scattering,  instead  of  the  highly  aniso¬ 
tropic  scattering  usually  found  in  nature,  or  the  dismissal  of  scat¬ 
tering  orders  greater  than  the  first  as  negligible. 

Coulson,  Dave  and  Sekera  (Reference  10)  published  a  rather  com¬ 
plete  set  of  calculated  data  on  the  transmission  and  reflection  of 
a  broad  beam  of  light  incident  on  the  top  of  a  planetary  atmosphere 
with  Rayleigh  scattering.  Their  method  of  calculation  was  that  de¬ 
rived  by  S.  Chandrasekhar  (Reference  11).  Although  the  calculational 
method  is  rigorous  for  the  model  considered,  it  is  not  applicable  to 
atmospheres  involving  scattering  processes  other  than  Rayleigh  scat¬ 
tering,  or  point  sources  of  light » 

It  appears  that  there  is  a  need  for  a  calculational  method 
which  does  not  require  the  use  of  fidelity-reducing  procedures  and 
which  would  provide  a  means  of  studying  the  effects  of  parameters 
such  as  cloud  cover,  ground  albedo,  and  atmospheric  aerosol  content 
on  the  transport  of  light  in  a  controlled  atmosphere.  Such  a  cal¬ 
culational  method  should  also  be  capable  of  studying  the  transport 
in  air  of  light  emitted  by  both  point  and  broad  beam  sources.  One 
method  that  has  the  capability  of  satisfying  all  of  the  above  re¬ 
quirements  is  the  Monte  Carlo  method. 
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The  Monte  Carlo  method  Is  distinguished  by  its  experimental 
nature.  It  involves  the  use  of  statistical  sampling  techniques  to 
approximate  the  solution  of  a  mathematical  or  physical  problem. 

The  application  of  the  Monte  Carlo  method  to  the  problem  of  light 
transport  in  the  atmosphere  makes  use  of  the  fact  that  the  radiation 
transport  equation  can  be  reduced  to  a  probabilistic  model  of  parti¬ 
cle  behaviour.  The  solution  of  the  transport  equation  can  then  be 
realized  by  averaging  the  results  of  individual  particle  histories, 
where  the  history  of  each  individual  particle  exhibits  a  random  be¬ 
haviour  in  accordance  with  the  probabilistic  model. 

A  preliminary  study  on  the  application  of  Monte  Carlo  methods  to 
the  light  transport  problem  was  performed  using  Monte  Carlo  codes 
that  had  originally  been  developed  for  the  study  of  neutron  trans¬ 
port.  The  results  of  this  preliminary  study  are  described  in  Section 
II  of  this  report. 

As  the  result  of  the  preliminary  study  on  the  application  of 
Monte  Carlo  methods  to  light  transport  problems,  it  became  apparent 
that  the  Monte  Carlo  method  could  provide  a  means  of  treating  rigor¬ 
ously  the  multiple  scattering  of  light  by  both  molecules  and  aerosol 
particles  in  a  slab  geometry  atmosphere.  Two  machine  programs  were 
developed  which  have  been  designated  as  LITE-I  and  LITE-I I.  LITE-I 
allows  for  treatment  of  problems  involving  point  sources  and  LITE-II 
treats  problems  involving  broad  beam  sources.  The  LITE  codes  were 
originally  written  in  FORTRAN-II  language  for  the  IBM-7090  computer. 
Following  checkout  and  debugging  of  the  FORTRAN  codes,  they  were 
converted  to  ALGOL  language  codes  by  the  use  of  a  FORTRAN-to-ALGOL 
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Translator  for  the  Burroughs  B-5000  computer.  In  Section  III,  a  de¬ 
tailed  presentation  of  the  methods  utilized  In  the  LITE  codes  is 
given.  Section  IV  shows  comparisons  of  LITE  codes  results  with 
several  other  sets  of  calculated  data. 

A  listing  of  the  ALGOL  language  codes  and  the  utilization  In¬ 
structions  for  the  codes  are  presented  In  Volume  II  of  this  report. 
Also  presented  In  Volume  II  are  listings  of  a  sample  problem  Input 
and  output  for  each  of  the  LITE  codes. 
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II.  PRELIMINARY  INVESTIGATIONS 

Th«  objective  of  the  preliminary  investigations  was  to  show  that 
the  transport  equation  describing  light  transport  in  the  atmosphere 
is  amendable  to  solution  by  Monte  Carlo  techniques  through  the  use  of 
existing  Monte  Carlo  codes  which  were  developed  for  the  purpose  of 

i 

studying  neutron  and  gamma-ray  transport.  The  results  of  the  Monte 
Carlo  calculations  obtained  frc.i  the  use  of  the  nuclear  radiation 
transport  codes  were  to  provide  information  on  the  important  para¬ 
meters  to  be  considered  in  the  formulation  of  a  Monte  Carlo  code 
which  was  to  be  developed  for  the  purpose  of  studying  light  trans¬ 
port  in  an  air-ground  geometry. 

An  initial  study  was  conducted  to  investigate  the  effect  of  higher 
order  scattering  on  the  total  scattered  intensity  for  a  point  source 
and  receiver  located  in  an  air-ground  geometry.  The  study  was  con¬ 
ducted  using  the  L05  Monte  Carlo  code  (Reference  13)  which  allowed  a 
description  of  an  air-ground  geometry  by  considering  the  ground  as  a 
dense  layer  of  atmosphere  in  which  scattering  was  assumed  isotropic 
in  the  laboratory  system. 

The  change  in  atmospheric  density  was  approximated  with  ten  slab 
layers  of  atmosphere  with  decreasing  densities.  The  extinction  cross 
section  for  air  at  ground  level  was  taken  to  be  4.403  x  10-^  cm-*. 

The  Rayleigh  and  Mie  scattering  cross  sections  were  taken  to  be 
2.523  x  10  ^  cm  *  and  1.88  x  10  ^  cm  *  respectively.  The  phase  func¬ 
tion  for  Mie  scattering  was  taken  to  be  that  computed  for  an  aerosol 

_3 

particle  size  distribution  proportional  to  r  ,  where  r  is  the  parti¬ 
cle  radius.  The  extinction ,  Rayleigh  scattering,  and  Mie  scattering 
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cross  sections  were  allowed  to  vary  in  the  same  manner  as  the  atmos¬ 
pheric  density  varys  with  altitude.  The  L05  code  was  used  to  com¬ 
pute  the  scattered  light  Intensity  from  a  point  isotropic  light 
source  located  1  km  above  the  ground  surface  for  a  point  receiver 
located  at  ground  level  10  km  horizontal  dist  nee  from  the  source 
(a  source-receiver  separation  distance  of  4.4  mean-free-path  lengths). 

A  ground  albedo  of  0.5  was  approximated  in  the  calculations  by  select¬ 
ing  the  proper  scattering-to-total  cross  section  ratio  for  the  ground 
layer.  The  scattering-to-total  cross  section  ratio  used  in  the  L05 
calculations  for  the  ground  layer  was  obtained  from  calculations  re¬ 
ported  by  Rafalski  ( Reference  14)  which  gives  the  neutron  albedo  as 
a  function  of  the  scattering-to-total  cross  section  ratio  for  a 
medium  in  which  scattering  is  isotropic  in  the  laboratory  system  and 
in  which  no  energy  degradation  is  allowed.  A  plot  of  Rafalski' s 
calculations  of  the  albedo  versus  the  scattering-to-total  cross  section 
ratio  is  shown  in  Figure  1.  Figure  2  shews  the  results  obtained  from 
five  L05  problems  in  which  the  maximum  number  of  collisions  allowed 
per  history  was  taken  as  5,  10,  20,  30  and  40,  respectively.  The 
results  obtained  from  problems  with  sample  sizes  of  1000  histories 
show  that  the  scattered  intensity  at  the  receiver  position  increases 
with  an  increase  in  the  maximum  number  of  collisions  allowed  per 
history  up  through  a  maximum  of  approximately  20  collisions  per  his¬ 
tory.  The  scattered  intensities  for  problems  involving  more  than  20 
collisions  per  history  are  approximately  equal.  The  calculated  point 
shown  in  Figure  2  for  a  maximum  of  20  collisions  per  history  has  a 
large  standard  deviation  and  is  felt  to  be  too  high,  particularly 
since  the  scattered  intensities  for  30  and  40  maximum  collisions  per 
history  are  both  below  the  value  computed  for  20  collisions. 
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Fig.  1.  Dependence  of  the  Scattering-to-Total  Cross  Section  Ratio  on  the  Albedo 
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Some  additional  information  on  the  five  L05  problems  is  given 
in  Table  I.  It  is  noted  that  although  the  average  number  of  colli¬ 
sions  per  history  increases  as  the  maximum  number  allowed  increases, 
the  average  number  of  collisions  per  history  appears  to  be  reaching 
a  maximum  value  as  the  number  of  maximum  collisions  allowed  per  his¬ 
tory  increases  to  40.  The  fact  that  a  maximum  value  in  the  «  verage 
number  of  collisions  followed  per  history  will  be  reached  is  indicated 
by  the  increase  in  the  number  of  histories  terminated  by  particle 
escape  from  the  problem  geometry  (leakage  into  the  atmosphere  above 
30  km).  From  the  results  shewn  in  Figure  2  and  Table  I,  it  appears 
that  valid  results  for  scattered  intensities  can  be  obtained  from 
L05  calculations  in  which  the  maximum  number  of  collisions  allowed 
per  history  is  20,  provided  that  an  adequate  number  of  histories  are 
treated. 

Since  the  maximum  number  of  collisions  that  should  be  considered 
in  any  given  problem  is  a  function  of  the  geometry  of  the  particular 
problem,  it  was  decided  that  the  light  transport  code  to  be  developed 
should  provide  information  on  the  scattered  intensity  resulting  from 
each  order  of  scattering.  This  will  allow  an  inspection  of  the  re¬ 
sults  to  determine  whether  the  intensity  for  the  highest  order  of 
collision  allowed  adds  significantly  to  the  total  scattered  intensity. 

In  order  to  determine  the  geometric  effects  of  an  air-ground  geo¬ 
metry  on  the  transport  of  light,  two  additional  L05  problems  were  run 
for  a  point  isotropic  light  source  where  the  source-receiver  separa¬ 
tion  distances  were  1.5  and  3  mean-free-path  lengths  respectively. 

The  results  of  the  three  L05  problems  were  then  compared  ;.*ith  the 
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Table  I.  Dependence  of  Various  Parameters  on  Maximum  Number  of 
Collisions  Allowed  Per  History  in  L05  Problems 


Maximum  No. 
of  Collisions 
Allowed  per 
History 

Fraction  of 
Histories 
Terminated 
by  Geometry 

Fraction  of 
Histories 
Terminated 
by  Max.  Coll. 

Average  No. 
of  Collisions 
per  History 

Intensity 
of  Scattered 
Light 

5 

0.474 

0.526 

3.239 

6.028-7  * 

10 

0.623 

0.377 

5.649 

1.785-6 

20 

0.707 

0.293 

9.312 

2.891-6 

30 

0.769 

0.231 

11.788 

2.130-6 

40 

0.822 

0.178 

12.968 

2.094-6 

*  Read  6.028-7  as  6.028  x  10 
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total  Intensities  calculated  for  the  respective  source-receiver 

separation  distances  for  the  source  and  receiver  embedded  within  an 

infinite  media.  The  Infinite  media  data  were  obtained  using  the  K74 

Monte  Carlo  code  (Reference  15)  and  the  properties  of  the  infinite 

media  were  identical  to  those  of  the  lower  layer  of  the  atmosphere 

in  the  L05  problems.  Figure  3  shows  a  comparison  of  the  L05  and  the 

2 

K74  data.  The  intensities  shown  have  been  multiplied  by  4nR  ,  where 

R  is  the  slant  range  between  the  source  and  receiver.  The  differences 

between  the  K74  and  the  L05  results  should  be  a  measure  of  the  effect 

of  a  varying  atmospheric  density  and  of  the  air-ground  interface  on 

2 

light  transport  in  air.  The  direct  intensity  times  4ttR  is  also  shown 
in  Figure  3  and  was  determined  to  be  the  same  with  both  codes.  The 
differences  between  tne  K74  and  L05  total  intensity  curves  and  the 
direct  intensity  curve  are  the  scattered  intensities  received  by  an 
isotropic  receiver. 

Since  the  geometric  effects  of  an  air-ground  geometry  greatly 
affect  the  scattered  light  intensity,  it  was  decided  that  the  light 
transport  code  to  be  developed  should  print  the  scattered  intensity 
in  a  form  that  would  allow  the  geometric  effects  to  be  analyzed.  In 
particular,  the  scattered  Intensities  should  be  printed  as  a  function 
of  both  the  region  of  scatter  and  the  number  of  times  a  particle  has 
reflected  from  a  given  reflection  surface. 

The  importance  of  multiple  scattering  was  also  investigated  for 
a  broad  beam  source  of  light  incident  upon  an  infinite  plane-parallel 
Rayleigh  atmosphere  of  5  mean-free-paths  in  thickness.  This  study 
was  conducted  using  the  C-18  Monte  Carlo  slab  penetration  code 


4nR7  TIMES  FLUX  (4ircm2  photons  cm“2/$ource  photon) 


13 


Fig.  3. 


Comparison  af  Calculations  af  Photon  Flux  Times  Distance  Squared  for  Infinite  Medium  and  Air-Graund 
Geometries 
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(Reference  16).  Light  from  a  broad  beam  source  was  taken  to  be  inci¬ 
dent  normal  to  the  atmosphere,  and  the  total  intensity  was  calculated 
at  various  depths  in  the  slab.  Figure  4  presents  the  results  obtained 
for  the  total  intensity  after  a  given  number  of  collisions  as  a  func¬ 
tion  of  depth  into  the  atmosphere.  The  importance  of  multiple  scat¬ 
tering  is  seen  to  increase  as  a  function  of  depth  into  the  atmosphere. 
The  intensity  after  30  collisions  is  only  a  factor  of  1.6  greater  than 
the  intensity  after  5  collisions  for  a  depth  of  1  mean-free-path, 
whereas  the  intensity  after  30  collisions  is  a  factor  of  4  greater 
than  that  after  5  collisions  for  a  depth  of  5  mean-free-paths. 

A  set  of  three  L05  problems  were  run  to  compare  with  light  scat¬ 
tering  data  supplied  by  USAEL  which  were  taken  during  an  Everglades, 
Florida  experiment.  Although  detailed  information  concerning  the 
ground  albedos  for  this  data  were  not  available,  the  albedos  were 
estimated  in  order  that  a  rough  comparison  might  be  made.  The  three 
L05  problems  were  designed  to  calculate  the  scattered  and  direct  in¬ 
tensities  received  by  a  flat  plate  receiver  at  distances  of  1,  2,  3, 

5  and  7  miles  for  both  blue  and  green  light.  The  light  source  was 
located  70  feet  above  ground  level  and  the  receiver  points  were  located 
12  feet  above  the  ground.  The  extinction  coefficients  used  for  each 
wave  length,  as  well  as  the  phase  functions,  were  obtained  from  the 
measured  data.  The  scattering-to-total  cross  section  ratio  of  the 
ground  media  was  selected  to  simulate  a  ground  albedo  of  approximately 
0.5.  The  results  for  the  green  light  source  are  presented  in  Figure 
5,  where  they  are  compared  with  measured  data.  Since  the  source  in¬ 
tensity  was  not  given,  the  magnitude  of  the  measured  intensities  were 


TOTAL  FLUX  (photons  cm“2/s°urc*  photon  cm 
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Fig.  4.  Variation  of  the  Total  Flux  After  a  Given  Number  of  Collisions  with  Penetrotion  Distance 


PHOTON  FLUX  (photons  cm' 
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RANGE  (miles) 


Fig.  5.  Light  Intensity  Received  by  a  Flat  Plate  Receiver  as  a  Function  af  Slant  Range:  Ground  Albedo  =  0.5 
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normalized  to  the  calculated  direct  Intensity  at  the  3  mile  distance. 
The  agreement  between  the  measured  and  calculated  data  are  rather  good, 
except  for  the  point  at  seven  miles  where  the  calculated  data  greatly 
overpredicts  the  measured  intensities.  It  is  thought  that  the  ground 
albedo  of  0.5  is  probably  too  high  for  the  Everglades  terrain.  A  high 
albedo  would  have  a  greater  effect  on  the  7  mile  point  than  on  the 
other  points.  Figure  6  shows  a  comparison  of  the  calculated  and  mea¬ 
sured  data  for  blue  light.  The  results  of  an  additional  calculation 
for  a  ground  albedo  of  zero  are  also  shown  in  Figure  6.  It  is  seen 
that  the  data  for  the  zero  ground  albedo  agree  better  with  the  mea¬ 
sured  intensities  than  do  the  results  for  the  0.5  ground  albedo.  The 
overall  agreement  between  the  two  sets  of  measured  and  calculated 
data  lends  confidence  in  the  use  of  Monte  Carlo  methods  to  predict 
scattered  intensity  of  light  in  an  air-ground  geometry. 


PHOTON  FLUX  (photons  cm" 
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Fig.  6.  Comparison  of  Calculations  of  Light  Transport  in  Which  the  Ground  Albedo  was  Taken  to  be  0,5  and  0.0 
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III.  METHODS 

The  favorable  results  obtained  with  the  nuclear  radiation  shield¬ 
ing  codes  L05,  K74  and  C18  during  the  preliminary  Investigations  of 
light  transport  in  the  atmosphere  led  to  the  development  of  two  Monte 
Carlo  procedures  designed  specifically  ter  the  study  of  light  trans¬ 
port.  The  first  of  these  procedures  was  developed  for  use  in  studying 
the  transport  of  light  emitted  isotropically  or  with  an  arbitrary 
polar  angular  distribution  by  a  point  source  located  in  an  air-ground 
geometry.  This  procedure  has  been  designated  as  the  LITE-I  code.  The 
second  procedure  (LITE-I I)  was  developed  to  study  the  transport  of 
light  emitted  from  a  plane  source  with  an  arbitrary  angular  distribu¬ 
tion  located  either  at  the  top  of  an  air-ground  atmosphere  or  within 
the  atmosphere. 

The  codes  treat  multiple  scattering  of  monochromatic  light  in  an 
atmosphere  that  may  vary  arbitrarily  with  altitude  Routines  are  avail¬ 
able  for  treating  both  Rayleigh  and  Mie  scattering  events.  An  inter¬ 
mixture  of  the  two  events  is  possible  or  an  atmosphere  may  be  con¬ 
sidered  as  either  a  Rayleigh  or  a  Mie  atmosphere.  The  Mie  scattering 
phase  functions  are  read  as  input  for  any  given  problem  so  the  aerosol 
concentrations  may  be  varied  from  one  problem  to  another.  The  LITE 
codes  do  not  consider  polarization  of  light 

Albedo  techniques  are  incorpoiated  to  treat  both  ground  and  cloud 
reflection;  however,  either  the  ground  or  cloud  regions  may  be  treated 
as  regions  of  increased  density  rather  than  with  the  albedo  method,  if 


desirable. 
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Basically ,  the  Monte  Carlo  codes  treat  light  on  a  particle  basis 
such  that  the  weight  assigned  to  each  individual  light  particle  re¬ 
presents  either  groups  of  photons  or  a  fractional  part  of  a  photon. 

The  path  that  an  individual  particle  traverses  in  an  air-ground  medium 
is  generated  by  a  random  selection  of  the  particle's  directions  and 
distances  between  collisions.  The  initial  direction  is  chosen  from 
the  input  source  angular  distribution  and  the  change  in  direction  at 
each  collision  is  determined  by  sampling  a  scattering  angle  from 
either  the  Rayleigh  or  Mie  distribution  depending  upon  which  scattering 
process  the  particle  has  undergone.  Distances  between  collisions  are 
chosen  by  selecting  a  path  length  from  either  an  exponential  or  a 
truncated  exponential  distribution  and  calculating  the  distance  corre¬ 
sponding  to  that  path  length.  A  weight  parameter  initialized  to  a 
value  of  one  is  assigned  to  each  source  particle.  Each  time  sampling 
from  a  biased  distribution  is  conducted,  the  weight  parameter  is  ad¬ 
justed  to  remove  the  bias  introduced. 

At  each  collision  an  estimate,  based  upon  the  particle's  weight, 
of  the  light  intensity  that  would  scatter  directly  to  each  of  a  set 
of  one  to  ten  receiver  positions  is  made.  Point  omnidirectional  re¬ 
ceivers  are  used  in  the  L1TE-I  code  and  plane  receivers  are  used  in 
the  LITE-II  code.  One  history,  composed  of  the  tracing  of  an  in¬ 
dividual  particle  through  all  orders  of  scattering  until  it  is  ter¬ 
minated,  provides  an  estimate  of  the  total  scattered  light  intensity 
at  each  of  the  receiver  positions.  A  final  estimate  of  the  scattered 
intensities  is  obtained  by  averaging  the  results  for  each  of  the  in¬ 


dividual  histories. 
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In  order  to  examine  the  statistical  fluctuation  of  the  results, 
the  total  number  of  histories  may  be  subdivided  into  groups,  each 
containing  an  equal  number  of  histories  such  tha*  any  one  group  does 
not  contain  data  for  more  than  500  histories  The  average  intensities 
for  each  individual  group  of  histories  are  printed  as  a  function  of 
the  order  of  collision  and  detector  numbers .  A  visual  examination  of 
the  statistical  fluctuation  is  then  possible  by  plotting  the  results 
for  each  of  the  groups. 

Upon  completion  of  the  last  group  of  histories  the  intensities 
at  each  receiver  as  a  function  of  order  of  collision  are  averaged 
over  the  total  number  of  groups  and  then  printed.  The  deviation  of 
the  group  results  about  this  average  is  also  calculated  and  printed. 
The  scattered  intensities  at  each  receiver  position  are  then  printed 
as  a  function  of  angle  and  order  of  reflection  for  a  given  reflection 
surface.  The  scattered  intensities  are  also  printed  as  a  function  of 
region  of  scatter  to  aid  in  the  determination  of  those  regions  which 
contribute  most  significantly  to  the  scattered  intensity. 

In  order  to  examine  the  amount  of  light  reflected  from  all  of 
the  reflection  surfaces  to  each  receiver  position,  the  reflected 
light  intensities  for  each  receiver  are  also  printed.  After  the  total 
scattered  intensities  are  calculated  and  printed,  the  direct  inten¬ 
sity  at  each  receiver  position  is  calculated  and  printed. 

The  following  portions  of  Section  £11  describe  in  detail  the 
methods  utilized  in  the  LITE  codes  The  formulae  presented  are  a 
mixture  of  the  usual  mathematical  notation  and  that  used  in  FORTRAN 
language.  The  *  is  used  to  denote  multiplication  and  parentheses 
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are  often  used  to  enclose  subscripts.  The  LITE-I  code  is  discussed 
initially  with  a  following  discussion  of  those  features  in  the  LITE-II 
code  that  differ  from  the  LITE-I  code. 


3.1  Description  of  LITE-I  Code 
3.1.1.  Selection  of  Source  Annies 


The  LITE-I  code  provides  for  the  description  of  a  point  isotropic 
or  point  anisotropic  source.  The  source  is  located  on  the  vertical 
axis  at  a  height  HS  above  the  air-ground  interface.  The  angular  dis¬ 
tribution  of  the  source  is  defined  in  terms  of  the  cosine  of  the  polar 
angle  measured  from  the  vertical  H  axis.  It  is  assumed  that  the  angu¬ 
lar  distribution  is  azimuthally  symmetric.  The  input  values  PAG(J) 
defining  the  polar  angular  distribution  are  given  by  the  expression 


PAG(J) 


CANG(J) 

-n(cos0)dcos0 

CANG(l) _ 

CANG(NAG) 

-n(cos0)dcos0 

CANG(l) 


t 


where  n(cos0)  is  the  function  describing  the  desired  angular  distri- 
but ion , 

CANG(J)  are  the  input  cosine  values  used  in  defining  the  angular 
distribution, 

NAG  is  the  number  of  input  cosine  values,  and 
PAG(l)  -  0  for  J  -  1. 


If  the  distribution  is  isotropic,  n(cos 0)  *  1  and 
PAG(J)  =  “  (CANG(l)  -  CANG(J) ) 
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Selection  of  Cosine  of  Source  Polar  Angle 


PAG(J)  -  PAG(J-l)  is  the  probability  of  the  particle  being 
emitted  within  the  interval  from  CANG(J-l)  to  CANG(J).  The  interval 
through  which  a  particle  is  emitted  is  selected  at  random  by  locating 
the  two  values  of  PAC(J)  such  that  PAG(J-l)  iRNl  <PA.G(J),  where  RN1  is 
a  random  number.  The  cosine  of  the  emission  angle  is  obtained  by  a 
linear  interpolation  with  a  second  random  number  between  the  two  values 
of  CANG(J)  that  bound  the  interval 

SANG  «  CANG(J-l)  -  RN2  (CANG(J-l)-CANG(J) ) 
where  SANG  is  the  cosine  of  the  source  angle  selected, 

RN2  is  a  second  random  number,  and 
CANG(J)  are  the  input  cosines  bounding  the  interval  selected  with 


the  first  random  number. 
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For  some  problems  It  may  be  desirable  to  bias  the  sampling  from 
the  source  polar  angular  distribution  to  improve  the  sampling  in  those 
directions  which  point  toward  the  detector  positions.  The  biased  dis¬ 
tribution  is  input  into  the  code  and  the  weight  of  the  light  particle 
is  adjusted  for  its  not  having  been  sampled  from  the  true  angular  dis¬ 
tribution.  If  the  true  angular  distribution  is  isotropic  (NAOP  *  0) , 
then  the  code  adjusts  the  particle  weight  internally  with  the  expres¬ 
sion 


WEIGHT (I) 


r _ LJ) _ -v /CANG(X-l)  -  PANG  (I)  ^ 

'■PAG (I)  -  PAG(I-1)^CANG(1)  -  CANG(NAG)  ; 


If  the  true  distribution  is  anisotropic  and  a  biased  distribution 
is  input  (NAOP  *  1),  the  weight  adjustment  factors  WAG(J)  must  also 
be  input.  WAG(I)  is  the  weight  that  is  assigned  to  all  particles 
emitted  from  the  source  with  the  cosines  of  their  angles  between 
CANG(I-l)  and  CANG(I).  WAG(l)  is  not  used  by  the  code  since  PAG(l) 
is  equal  to  0  and  the  code  never  selects  an  I  =  1.  WAG (2)  is  the 
weight  that:  will  be  assigned  to  particles  emitted  from  the  source 
with  a  cosine  between  CANG(l)  and  CANG(2).  The  WAG(J)  values  should 
be  computed  with  the  following  equation, 


WAG(J) 


CANG(J)  ,,  ...  u 

j  -n  (cos0)dcos0 

CANG(J-l) _ 

CANG(J) 

/  -n(cos0)dcos0 

CANG(J-l) 


where  n* (cos0)  is  the  true  angular  distribution,  and 
n(cos0)  is  the  input  biased  distribution. 


The  LITE-I  code  presently  treats  only  those  source  angular  dis¬ 
tributions  that  are  azimuthally  symmetric.  Rather  than  selecting  the 
azimuthal  angle  from  an  isotropic  distribution,  all  source  angles  are 
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started  with  an  azimuthal  angle  of  0°.  Each  time  an  estimate  of  the 
scattered  intensity  at  a  receiver  is  made,  the  azimuthal  position  of 
the  receivers  are  selected  from  a  uniform  distribution  between  0  and 
2n.  This  produces  the  same  results  as  selecting  a  source  azimuthal 
angle  from  an  isotropic  distribution  and  estimating  the  light  inten¬ 
sity  at  a  receiver  whose  azimuthal  position  is  fixed.  The  efficiency 
of  the  LITE-I  code  may  be  increased  by  selecting  two  or  more  azimuthal 
positions  for  a  receiver  point  and  averaging  the  contributions  to  all 
the  positions  from  each  collision.  The  process  of  selecting  N  re¬ 
ceiver  azimuthal  positions  and  averaging  the  contributions  to  each 
position  is  equivalent  to  increasing  the  original  number  of  histories 
by  a  factor  of  N.  This  process,  however,  results  in  a  considerable 
amount  of  machine  time  3aved  since  the  same  set  of  collision  points 
is  used  N  times. 

3.1.2  Optical  Path  Lengths 

Since  the  extinction  coefficient  for  a  clear  atmosphere  has  been 
shown  to  differ  considerably  from  the  exponential  variation  at  the 
upper  end  of  the  visible  range  and  also  below  the  visible  range  (Ref. 
17),  a  method  was  developed  to  allow  an  arbitrary  variation  of  extinc¬ 
tion  coefficient  with  altitude.  Rather  than  inputting  the  extinction 

coefficient,  the  optical  thickness  from  ground  level  to  altitude  hi, 
hi 

t (hi)  =  /  E(h)dh, 

0 

is  input  for  values  of  hi  sufficiently  close  enough  together  so  as  to 
adequately  represent  the  curve  t(h)  with  straight  line  segments  join¬ 
ing  the  i(hl)  values.  £(h)  is  the  extinction  coefficient  at  altitude  h. 
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For  a  light  particle  located  at  an  altitude  hi  and  traveling  in  a 
direction  making  an  angle  6  with  the  vertical  vector,  the  number  of 
path  lengths,  RHO,  traversed  before  collision  is  sampled  from  an  ex- 
ponential  distribution  if  the  particle's  direction  is  downward  or 
from  a  truncated  exponential  distribution  if  the  particle's  direction 
is  upward.  Sampling  from  the  exponential  distribution  is  accomplished 
by  solving  the  integral  equation 
RHO 

RN  -  J  e“pdp  , 

0 

for  RHO,  when  RN  is  a  random  number  between  0  and  1.  For  the  trun¬ 
cated  exponential  distribution  the  path  length,  RHO,  is  given  by  the 
equation, 

RHO 

RN  =  /  e~pdo  3. 1.2.1 

®  l-e-pmax 

where  Pmav  is  the  number  of  path  lengths  from  the  collision  point  at 
hi  to  the  upper  surface  of  the  atmosphere  along  the  particle's  direc¬ 
tional  vector.  The  quantity  Pmax  is  given  by  the  expression 

p  =  (x(h  )  -  x(hl))  /co80  3.1. 2.2 

where  x (hmay)  is  the  thickness  of  the  atmosphere  in  mean-free-path 
lengths, 

x(hl)  is  the  number  of  mean-free-path  lengths  from  the  ground 
to  altitude  hi,  and 

cos0  is  the  cosine  of  the  angle  that  the  particle's  direc¬ 
tional  vector  makes  with  the  positive  vertical  vector. 
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The  solution  of  equation  3. 1.2. 2  for  RHO  gives 

RHO  -  -  f.n(l«0  -  RN(1.0 -e"0®**)} 

which  reduces  for  p  •  •  to 

max 

RHO  -  -  £n(1.0  -  RH)  , 

the  solution  of  equation  3. 1.2.1,  the  case  when  the  particle's  direc¬ 
tion  is  downward. 

Sampling  path  lengths  from  the  truncated  exponential  distribu¬ 
tion  introduces  biasing  since  a  collision  is  forced  before  the  particle 
escapes  the  upper  atmosphere.  This  biasing  is  removed  by  multiplying 
the  particle's  weight  by 
1  -  e~Pmax, 

the  probability  a  collision  will  occur  before  it  escapes  the  atmosphere. 

Distance  Between  Collision 


h 

max 


DTST*sin0 
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After  the  path  length  RHO  has  been  sampled,  the  altitude  of  the 
next  collision  Is  determined  by  solving  the  expression 

t(h2)  »  t (hi)  +  RHO*cos0 

and  Interpolating  between  the  values  of  hi  corresponding  to  the  two 
values  of  t(hl)  that  bound  t(h2)  to  find  the  height  h2  that  corre¬ 
sponds  to  t(h2). 

The  distance  between  collisions  Is  computed  with  the  formula 
DIST  *  (h2  -  hi) /cos0 

where  hi  Is  the  altitude  of  nth  collision, 

h2  Is  the  altitude  of  the  (n+l)th  collision,  and 
cos 6  Is  the  cosine  of  the  angle  between  the  line  joining  the  two 
collision  points  and  the  positive  vertical  vector. 

The  radial  position  of  the  (r.+l)th  collision  point  is  computed 
with  the  expression 

R2  =t/r12  +  (DIST*sine)2  +  2*Rl*DIST*sin6*cos$ 

where  Rl  is  the  radial  position  of  the  nth  collision, 

DIST*sin0  is  the  projection  of  the  line  joining  the  nth  and  (n-H)th 
collisions  into  the  horizontal  plane,  and 
cos$  is  the  angle  that  the  horizontal  projection  makes  with  a 
radial  vector  through  the  nth  collision  point. 

3.1.3  Atmosphere  Description 

Since  the  aerosol  size  distribution  and  the  variation  of  aerosol 
concentration  with  altitude  may  vary  independently  with  air  density, 
the  ratios  of  the  Rayleigh- to-scattering  cross  section  and  the  ratios 
of  the  scattering -to-total  cross  section  are  input  for  the  points  hi 
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used  in  defining  t(h).  The  probability  of  Mie  scattering  at  any 
altitude  is  then  given  by 
1  -  RAYR(h) 

where  RAYR(h)  is  the  ratio  of  the  Rayleigh-to-scattering  cross  section. 

The  variation  of  aerosol  particle  size  distribution  with  alti¬ 
tude  is  treated  by  dividing  the  atmosphere  into  regions  and  defining 
a  different  Mie  scattering  phase  function  for  each  region.  The  Mie 
scattering  angle  for  each  Mie  scattering  event  is  chosen  from  the 
phase  function  for  the  region  containing  the  collision  point. 

Light  scattering  within  cloud  layers  can  be  considered  by  the 
code,  if  desired.  If  a  cloud  layer  is  located  between  heights  HC1 
and  HC2,  where  HC1  and  HC2  are  the  heights  of  the  cop  and  bottom  of 
the  cloud  layer,  then  the  extinction  coefficient  E(h)  used  in  the 
equation 

hi 

t (hi)  -  /  E(h)dh 
0 

for  HCl>hiHC2  should  be  that  for  the  cloud  layer.  The  values  of 
the  ratio  RAYR(h)  input  for  HCl.2h2.HC2  should  be  the  ratio  of  the 
Rayleigh  to  the  total  scattering  cross  section  for  the  cloud  layer. 

The  cloud  layer  can  be  defined  by  one  or  more  regions  and  different 
Mie  scattering  phase  functions  can  be  input  for  each  region,. 

3.1.4  Geometry 

A  description  of  an  air-ground  geometry  is  accomplished  by 
dividing  the  geometry  into  regions  bounded  by  horizontal  plane  sur¬ 
faces 


0  =  COEE(NCB)  -  H 


3. 1.4.1 
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and/or  right  circular  vertical  cylindrical  surfaces 

0  =  COLj(NCB)  -  R  ,  3. 1.4. 2 

where  H  and  R  are  coordinates  of  a  point  on  one  of  the  surfaces, 
COEE(NCB)  is  the  H  intercept  of  the  plane  surface  or  the  radius  R 
of  the  cylindrical  surface,  and  NCB  is  boundary  number  for  a  given 
surface.  A  region  is  defined  by  listing  the  numbers  of  the  boundary 
surfaces  encompassing  the  region.  A  plus  or  minus  sign  is  affixed  to 
each  of  the  boundaries  listed  for  a  region  to  designate  the  boundary 
as  an  "upper"  or  "lower"  plane  surface  or  an  "outer"  or  "inner"  cy¬ 
lindrical  surface  with  respect  to  the  region  being  defined. 

3. 1.4.1  SEARCH  Routine 

In  locating  the  region  containing  a  particle's  position  H,R, 
the  SEARCH  routine  calculates  a  value  XR  for  each  boundary  encompass¬ 
ing  the  region  by  use  of  one  of  the  equations 
XR  =  COEE(NCB)  -  H 
XR  =  COEE(NCB)  -  R 

and  compares  the  sign  of  XR  with  the  sign  affixed  to  the  correspond¬ 
ing  boundary  number  in  the  region  input  data.  If  the  signs  agree 
for  all  boundaries  surrounding  the  region,  then  the  point  H,R  is  in 
the  region  under  consideration.  If  not,  then  the  SEARCH  routine 
chooses  the  next  region  listed  in  the  region  table  and  repeats  the 
process  for  that  region  The  SEARCH  routine  continues  until  all  re¬ 
gions  have  been  searched  or  until  the  region  containing  the  point  H,R 
has  been  located.  The  search  is  always  begun  with  the  most  probable 
region  of  entry,  MPR,  for  the  last  boundary  that  has  been  crossed. 

An  error  in  listing  the  most  probable  region  of  entry  will  not  re¬ 
sult  in  an  error  in  the  calculations,  but  will  Increase  the  running 
time  for  a  given  problem. 
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3. 1.4. 2  Distance  to  a  Boundary 

The  distance  along  the  line  of  the  particle’s  direction  to  a 
surface  bounding  the  region  is  calculated  in  subroutine  DSTBD.  Sub¬ 
routine  DSTBD  actually  calculates  the  distances  to  all  boundaries 
surrounding  the  region  and  takes  the  shortest  positive  distance  as  the 
distance  traversed  within  the  region.  At  any  point  a  particle's  posi¬ 
tion  is  given  by  the  coordinates  H,R  and  the  direction  is  given  by 
the  two  angles  0  and  <f>.  6  is  the  angle  the  particle's  direction 

makes  with  the  positive  vertical  vector  and  $  is  the  angle  the  projec¬ 
tion  of  the  particle's  direction  into  a  horizontal  plane  makes  with 
the  radial  vector  through  the  point  H,R. 

Distance  to  Plane  Surface 
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The  distance  X  to  a  plane  surface  is  calculated  with  the  expression 
X  =  (COEE(NCB)  -  H)/COTH, 

where  COEE(NCB)  is  the  H  intercept  of  the  plane  surface  number  NCB,  and 
COTH  *  cos0  is  the  cosine  of  the  angle  between  the  parti¬ 
cle's  direction  and  the  vertical  vector. 

The  distance  X  to  a  cylindrical  surface  is  given  by  the  expression 

x  »  -r*cphi  t  /coee2(ncb)  -  R2*SPHT^ 

SITH 

where  COEE(NCB)  is  the  radius  of  the  cylindrical  surface  NCB, 

CPHI  and  SPHI  are  the  cosine  and  sine  of  the  angle  that  the  horizon¬ 
tal  projection  of  the  particle's  direction  makes  with 
the  radial  vector  through  H,R,  and 


SITH  is  the  sine  of  the  angle  the  particle's  direction 
makes  with  the  vertical  vector. 

Distance  To  Cylindrical  Surface 
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If  the  value  under  the  radical  is  negative,  then  there  is  no  intersec¬ 
tion  with  the  cylindrical  surface  NCB  When  the  value  under  the  radi¬ 
cal  is  positive,  the  line  intersects  the  cylindrical  surfaces  in  two 
points  If  the  radical  position  cl  the  point  R  is  less  than  COEE(NCB), 
then  the  point  H,R  is  inside  the  cylindrical  surface  and  the  positive 
radical  is  used  to  compute  the  distance  X,  but  if  the  radial  position 
of  the  poin*-  is  greater  than  COEE(NCB),  then  the  point  H,R  is  outside 
the  cylindrical  surface  and  the  negative  radical  is  used  to  compute 
the  distance  X, 

The  distance  DIST  that  a  particle  may  travel  within  a  region  is 
initially  set  equal  to  the  input  value  DLONG.  Each  time  a  distance  X 
to  a  boundary  surrounding  the  region  is  calculated,  X  is  tested  to 
determine  whether  it  is  a  positive  oc  negative  quantity  If  X  is  nega¬ 
tive,  then  DIST  remains  unchanged,  but  when  X  is  positive,  a  compari¬ 
son  of  X  and  DIST  is  made  If  DIST  is  iess  than  X,  then  the  value  for 
DIST  remains  unchanged,  but  when  DIST  is  greater  than  X,  the  value  of 
DIST  is  set  equal  to  X  After  calculating  the  distances  to  all  bounda¬ 
ries  surrounding  a  region  and  performing  the  above  tests,  the  final 
value  of  DIST  should  be  the  distance  to  the  first  boundary  cf  intercept, 

3, 1,4, 3  Region  Importance  Numbers 

The  region  description  provided  by  the  LITE-1  code  not  only 
allows  for  the  input  of  different  Mie  scattering  phase  functions  for 
each  region  of  the  atmosphere  but  it  also  allows  the  user  to  emphasize 
the  sampling  of  collision  points  in  those  regions  felt  to  be  of  most 
importance.  Each  region  is  assigned  an  importance  number,  EMP,  in  the 
problem  input  data  As  a  Light  particle  passes  from  one  region  to 


34 


another*  the  importance  numbers  of  the  two  regions  are  compared.  If 
the  particle  enters  a  region  of  the  same  or  more  importance,  tracing 
of  the  particle  is  continued  in  a  normal  manner.  If  the  particle  enters 
a  region  of  less  importance,  the  ratio  of  the  importance  number  of  the 
region  entered  to  the  importance  number  of  the  previous  region  is  com¬ 
pared  with  a  random  number.  If  the  ratio  is  less  than  the  random  number, 
the  particle  is  terminated.  If  the  ratio  is  greater  than  the  random 
number,  the  particle  weight  is  divided  by  the  ratio  and  tracing  of  the 
particle  is  continued.  This  method  makes  it  possible  to  concentrate 
the  sampling  of  particle  collisions  in  those  regions  that  contribute 
most  significantly  to  the  scattered  intensity  at  the  receiver  positions. 
If  a  region  is  assigned  an  importance  number  of  zero,  then  the  particle 
histories  for  all  particles  entering  that  region  will  be  terminated. 

3.1.5  Statistical  Estimation 

The  LITE-I  code  uses  Monte  Carlo  methods  to  generate  a  random 
walk  that  gives  the  location  of  a  particle's  collisions  within  an 
air-ground  geometry.  A  statistical  estimation  process  is  applied 
to  predict  the  light  intensity  that  scatters  from  each  collision 
point  of  a  random  walk  to  each  of  a  set  of  from  one  to  ten  point  de¬ 
tectors.  Since  only  the  height  and  radius  of  each  collision  point 
are  preserved  to  identify  its  location,  the  azimuthal  position  <J>^  of 
the  detector  point  is  allowed  to  vary  with  each  individual  estimate 
to  account  for  the  three  dimensional  collision  distribution.  The 
code  presently  selects  the  azimuthal  position  of  the  detector 
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point  from  a  uniform  distribution  between  0  and  2",  which,  in  effect, 
restricts  the  source  angular  distribution  to  one  that  is  symmetrical 
about  the  vertical  axis.  The  collisions  are  all  assumed  to  lie  on 
the  <f>jj  ■  0  plane  and  the  angle  ^  between  the  collision  plane  and  the 
radius  through  the  detector  point  is  chosen  with  a  rejection  technique. 

Geometry  For  Calculation  of  Flux  at  a  Detector  Point 
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3. 1.5.1  Detector  Angle  of  Scatter 

The  cosine  CSA  of  the  angle  between  the  particle's  direction 
before  collision  and  the  line  joining  the  collision  and  detector 
points  is  computed  by  taking  the  dot  product  of  the  two  unit  vec¬ 
tors  in  those  directions.  The  unit  vector  in  the  particle's  direc¬ 
tion  is 

I  -  (SITH*CPHI)x  +  (SITH*SPHI) j  4  (COTH)£ 
where  SITH  and  COTH  are  the  sine  and  cosine  of  the  angle  between 

the  particle's  direction  and  the  H  axis,  and 
SPHI  and  CPHI  are  the  sine  and  cosine  of  the  angle  between 
the  horizontal  projection  of  the  particle's 
direction  and  the  radial  vector  through  the 
point  of  collision. 

The  unit  vector  in  the  direction  of  the  line  joining  the  collision 
and  detector  points  is 

where  SPHID  and  CPHID  are  the  sine  and  cosine  of  the  detector's 
azimuthal  angle, 

Rp  and  Hp  are  the  radius  and  height  of  the  detector 
point , 

R  and  H  are  the  radius  and  height  of  the  collision 
point. 

The  distance  T  from  the  collision  point  to  the  detector  point  is 
computed  by  use  of  the  expression 

T  =  n/r2  4  ^  -  2RRdcos4>d  4  (H  -  Hp)^  . 

CSA  is  then  determined  by  the  dot  product 
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3. 1.3.2  Probability  of  Scattering  Toward  a  Detector 

CSA  la  the  cosine  of  the  angle  through  which  a  light  particle 
must  scatter  In  order  to  be  directed  toward  the  detector.  Subrou¬ 
tine  DIFSCA  computes  the  probability  of  scattering  within  a  unit 
solid  angle  about  the  angle  whose  cosine  is  CSA  Two  different 
methods  are  used  to  compute  the  Rayleigh  and  Mie  scattering  proba¬ 
bilities.  To  determine  the  mode  of  scattering  a  random  number  is 
generated  and  compared  to  the  value  input  for  RAYR(Hl) ,  which  is 
the  ratio  of  the  Rayleigh  to  the  total  scattering  cross  section  at 
the  altitude  HI  of  the  collision  point.  When  the  random  number  is 
less  than  the  value  for  RAYR(Hl)  (Rayleigh  scattering  occurs) ,  the 
probability  of  scattering  through  an  angle  wnose  cosine  is  CSA  is 
given  by  the  expression 

PSCAT=-r|“  (10  +  CSA2)  . 

When  the  random  number  is  greater  than  the  value  of  RAYR(Hl) ,  the 
probability  PSCAT  is  chosen  from  the  input  Mie  scattering  phase 
function  for  the  region  NCM  containing  the  collision  point.  The 
interval  containing  the  cosine  value  CSA  is  located  and  a  linear 
interpolation  between  the  two  input  values  of  PDCOS(NCM,J)  corre¬ 
sponding  to  the  two  cosine  values  DIFCOS(NCM,J)  that  bound  the 
interval  is  performed  to  obtain  the  value  PSCAT  for  the  cosine  CSA. 

31-5,3  Contribution  From  A  S ingig  Collision 

The  attenuation  of  the  light  atteied  from  a  collision  point 
to  a  detector  point  is  determined  by  calculating  the  number  of  path 
lengths  RHOT  traversed  along  the  line  joining  the  collision  and 


38 


detector  point.  RHOT  is  computed  by  use  of  the  equation 
RHOT  -  (TAUHD  -  TAUH2)cos8 

where  TAUH2  is  the  collision  altitude  in  mean-free-path  lengths, 

TAUHD  is  the  detector  altitude  in  mean-free-path  lengths,  and 
cosO  is  the  cosine  of  the  angle  between  the  line  joining  the 
collision  and  detector  point  and  the  positive  vertical 
vector. 

The  contribution  to  the  light  intensity  from  a  single  collision  is 
given  by  the  expression 

RESULT  =  (WAIT*PSCAT*e_RH0T) /T2 
where  WAIT  is  the  particle  weight  before  scattering, 

PSCAT  is  the  probability  of  scattering  into  a  unit  solid  angle 
about  the  direction  toward  the  detector, 

RHOT  is  the  number  of  path  lengths  between  the  collision  and 
detector  points,  and 

T  is  the  distance  from  the  collision  to  the  detector  point. 
3.1.6  Light  Scattering 

Each  individual  scattering  event  is  assumed  to  be  either  Rayleigh 
or  Mie  scattering.  For  a  given  atmosphere  the  ratio  of  the  Rayleigh 
to  the  total  scattering  cross  section  RAYR(hl)  is  input  as  a  function 
of  altitude  where  the  hi  are  fixed  heights  in  the  atmosphere.  The 
scattering  process  is  selected  by  comparing  a  random  number  with 
RAYR(h),  where  RAYR(h)  is  found  by  interpolation  between  the  RAYR(hl) 
values.  When  the  random  number  is  less  than  RAYR(h),  the  routine  for 
Rayleigh  scattering  is  used  to  select  the  scattering  angle.  Otherwise 
the  Mie  scattering  routine  is  used. 
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3. 1.6.1  Rayleigh  Scattering 

Rejection  techniques  are  used  to  sample  the  polar  and  azimuthal 
scattering  angles  from  the  Rayleigh  phase  function  which  is  expressed 
by  the  formulae 

P(fl)  =  77”  (1  +  cos20)dn  3. 1.6.1 

ICTT 

where  dft  *  sin0d0d<f>. 


Equation  3. 1.6.1  xs  expressed  as  the  product  of  two  probability 
d<i  3  2 

density  functions,  ~  times  -g  (1  +  cos  0)sin0d0.  The  azimuthal  angle 
of  scatter  is  selected  from  the  probability  density  function 

p<*>  -  # 


and  the  polar  angle  0  is  selected  from  the  probability  density  function, 
P (cos0)  =  ■—  (1  +  cos20)sin0d0  . 


Selection  of  the  azimuthal  angle  is  accomplished  by  use  of  the 
following  well-known  rejection  technique.  An  X  is  chosen  from  a  uni¬ 
form  distribution  in  the  interval  -1<X<1  by  use  of  the  equation 
X  =  2*RN.-1, 

and  a  Y  is  chosen  from  the  same  distribution  by  use  of  the  equation 
Y  =  2*RN2-1, 

where  RN^  and  RN2  are  random  numbers.  If  the  point  X,  Y  lies  outside 

a  circle  with  unit  radius,  then  two  more  values  of  X  and  Y  are  chosen 

and  the  process  is  repeated  until  a  point  X,  Y  is  found  within  the 

circle.  The  sign  and  cosine  of  the  azimuthal  scattering  angle  are 

then  given  by  the  equations 
cos0  =  — £ - 

/x2+  Yi 

sin0  =  X 

/X2+  Y2. 
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A  rejection  technique  was  developed  to  select  a  polar  angle  6  from 

the  probability  density  function, 

3  2 

P(cos0)  «  g-  (1  +  cos  0)sin0d0  . 

For  this  technique  the  above  density  function  is  expressed  as  the  sum 
of  two  probability  density  functions 
P(cos0)  *  P^(cos0)  +  P2(cos0) 
where  P^CcosO)  *  ~ 

P2(cos0)  “  cos^0  si29d9  • 


The  density  function  P^(cos0)  is  chosen  to  represent  the  polar 
angular  distribution  (3/4) / (6/4)  or  of  the  time,  and  the  density  func¬ 
tion  P2(cos0)  is  chosen  the  remainder  of  the  time.  If  the  density 
function  P^(cos0)  is  chosen,  then  cosd  is  determined  by  solving  the 


integral. 


RN  •*  / 
0 


P  sinQdQ 
2 


for  cosQ  . 

P 

Thus  cos  6  »  1-2*RN  . 

P 

When  the  density  function  P2(cos6)  is  chosen,  cosQ^  is  again  computed 

by  the  above  expression,  but  the  value  of  cosS  is  not  accepted  unless 

2 

a  second  random  number  is  less  than  cos  3  .  If  cos3  is  not  accepted, 

P  P 

then  the  process  of  selecting  a  density  function  and  computing  a  value 

of  cos3  is  repeated  until  either  P, (cosQ)  has  been  used  to  obtain  a 
P  1 

value  of  cosQ  or  Po(cos0)  has  been  used  and  the  value  has  been  ac- 
P  l 

cepted.  The  expected  number  of  trials  before  accepting  a  cosQ^  is  3/2 
and  the  efficiancy  of  the  technique  is  2/3. 


A  FORTRAN  code  was  prepared  for  the  purpose  of  verifying  the 
accuracy  of  the  rejection  technique  described  above.  The  printout  of 
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the  code  gives  the  number  of  angles  sampled  within  each  .05  interval 
of  cos6  between  1  and  -1  Figure  7  shows  a  comparison  of  the  cumula¬ 
tive  distribution  of  the  scattering  angles  sampled  with  the  code  and 
the  cumulative  Rayleigh  phase  function.  The  good  agreement  between 
the  two  curves  indicates  the  rejection  technique  does  properly  sample 
scattering  angles  from  the  Rayleigh  phase  function. 

3. 1.6. 2  Mie  Scattering 

The  azimuthal  angle  for  the  Mie  scattering  process  is  chosen  in 
the  same  manner  as  that  described  for  Rayleigh  scattering.  The  Mie 
scattering  polar  angle  is  selected  from  the  input  cumulative  Mie 
scattering  angular  distribution.  For  each  region  J  there  are  NPHANG(J) 
cosine  values,  PHANG(J,I),  input  corresponding  to  equal  intervals  of 
the  accumulative  Mie  scattering  angular  distribution.  A  random  number 
is  generated  and  the  values  1-1  and  I  determined  such  that 
(I-1)/NPHANG(J)  iRN  £ I/NPHANG(J) . 

Cos0p  is  then  determined  by  a  linear  .interpolation  between  the  cosine 
values  PHANG(J,I-1)  and  PHANG(J,I)  to  give  the  cosine  of  the  angle 
correstn'-.-Ang  to  the  random  number. 

3. 1.6. 3  Direction  after  Collision 

The  light  particle's  direction  after  scatter  is  determined  from 
a  knowledge  of  its  direction  before  scattering  and  the  polar  and  azi¬ 
muthal  angles  of  scatter.  The  cosine  of  the  particle's  direction 
after  scatter,  C0TH1,  is  given  by  the  equation, 

C0TH1  -  COTH2*CSANG+SITH2*SSANG*CDPHI, 
where  C0TH1  is  the  cosine  of  the  angle  between  the  particle's  direc¬ 
tion  after  scatter  and  the  vertical  axis. 
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Comparison  of  Cumulative  Distribution  of  Angles  Sampled  Using  Rejection  Technique  with  Cumulative 
Rayleigh  Phase  Function 
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SITH2  and  C0TH2  are  the  sine  and  cosine  of  the  angle  between  the  parti¬ 
cle's  direction  before,  scattering  and  the  vertical  axis, 
SSANG  and  CSANG  are  the  sine  and  cosine  of  the  polar  angle  of  scatter, 
and  CDPHI  is  the  cosine  of  the  azimuthal  angle  of  scatter. 


Geometry  for  Calculation  of  Direction  After  Collision 
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The  sine  of  the  angle  between  the  particle's  direction  after  scatter¬ 
ing  and  the  vertical  axis  is  given  by 
SITH1  -  yi-O-COTHl2  . 

The  sine  and  cosine  of  the  angle  between  the  horizontal  projection  of 
the  particle's  direction  after  scattering  and  the  radial  vector 
through  the  collision  point  are  given  by  the  equations 
SPHIl  -  SPHI2*CDEPHI  +  CPHI2*SDEPHI ,  and 
CPHI1  -  CPHI2*CDEPHI-SPHI2*SDEPHI, 
where  SPHI2  and  CPHI2  are  the  sine  and  cosine  of  the  angle  between 

the  horizontal  projection  of  the  particle's  direction  be¬ 
fore  scattering  and  the  radial  vector  through  the  scatter¬ 
ing  point, 

SDEPHI  »  SSANG*SDPHI / SITH1 ,  and 

CDEPHI  *  (CSANG  -  C0TH1*C0TH2) /SITH1*SITH2  . 

3.1.7  Ground  and  Cloud  Reflections 

Light  reflection  from  the  ground  ar.d  cloud  surfaces  can  be  treated 
with  an  albedo  method.  If  the  albedo  method  is  to  be  used,  then  the 
reflection  surfaces  are  listed  first  in  the  boundary  table  of  the 
LITE -I  input  and  are  identified  by  the  negative  sign  preceding  the 
boundary  number  NBOUND.  The  value  assigned  to  the  parameter  ALBEDO(NCB) 
gives  the  fraction  of  the  incident  photon  current  that  is  reflected 
from  the  surface  NCB  and  the  value  input  for  the  parameter  JREFLT(NCB) 
determines  the  angular  distribution  of  the  albedo.  The  angular  dis¬ 
tribution  is  taken  to  be  isotropic  or  anisotropic  in  the  upper  hemi¬ 
sphere  if  JREFLT(NCB)  is  1  or  2,  or  isotropic  or  anisotropic  in  the 
lower  hemisphere  if  JREFLT(NCB)  is  3  or  4.  When  reflection  is  isotropic, 
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in  either  the  upper  or  lower  hemisphere,  no  input  reflection  distribu¬ 
tions  are  required.  If  reflection  is  anisotropic,  then  values  must  be 
input  to  define  both  the  reflection  distribution  and  the  cumulative  re¬ 
flection  distribution. 

3.1. 7.1  Reflection  Into  A  Detector 

The  code  treats  light  reflection  from  a  surface  as  if  the  point  of 
reflection  is  a  scattering  center  at  the  reflection  surface.  The 
probability  of  scattering  toward  the  detector  and  the  angle  of  scatter 
is  determined  from  the  reflection  distributions  rather  than  from 
Rayleigh  or  Mie  scattering  distributions.  The  input  reflection  dis¬ 
tributions  are  defined  about  the  normal  to  the  reflection  surface. 

Reflection  into  Detector 
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The  input  values  POR(NRfl,J)  define  an  anisotropic  reflection  distribu- 
uicr  <  ‘T  the  cosine  values  RFANG(NRB,J)  for  the  reflection  surface 
NRB.  Foi  anisotropic  reflection  the  probability,  PSCAT,  of  light  re¬ 
flecting  from  a  surface  into  a  unit  solid  angle  toward  a  detector  is 
determined  by  calculating  the  cosine,  CSA,  of  the  angle  between  the 
line  joining  the  scattering  center  to  the  detector  and  the  normal  to 
the  surface,  and  then  evaluating  the  probability  of  light  reflecting 
at  that  angle  by  a  linear  interpolation  between  the  two  values  of 
POR(NRE,J)  corresponding  to  the  values  of  RFANG(NRB,J)  that  bound  the 
interval  containing  CSA.  The  angular  distribution  of  the  current 
reflected  from  a  surface  is  assumed  to  be  dependent  only  on  the  polar 
angle  of  reflection  and  to  not  vary  with  the  azimuthal  angle  for  any 
given  polar  angle. 

For  isotropic  reflection  in  the  upper  or  lower  hemisphere  the 
probability,  PSCAT,  of  light  reflecting  into  a  unit  solid  angle  to¬ 
ward  the  detector  is  taken  by  the  code  to  be  1/2it. 

The  light  intensity  reflected  into  the  detector  from  the  scat¬ 
tering  center  is  given  by  the  expression, 

_purjT  2 

RESULT  =  WAIT*ALBEDO (NRB) *PSCAT*e  /T 
where  RESULT  is  the  contribution  to  the  light  intensity  from  the 
scattering  center  at  the  reflection  surface, 

WAIT  is  the  light  particle  weight  before  scattering, 
ALBEDO(NRB)  is  the  fraction  of  the  particle  that  is  reflected  from 
the  surface  NRB, 

PSCAT  is  the  probability  of  the  particle  being  reflected 
into  a  unit  solid  angle  toward  the  detector, 
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RHOT  is  the  number  of  optical  path  lengths  between  the  scatter¬ 
ing  center  and  the  detector,  and 
T  Is  the  distance  between  the  scattering  center  and  the  de¬ 
tector. 

3.1. 7.2  Reflection  Angle 

After  computing  the  light  intensity  reflected  to  the  detector, 
a  reflection  angle  is  chosen  from  the  cumulative  reflection  distribu¬ 
tion  in  order  to  trace  the  light  particle  to  the  next  collision 
point.  If  reflection  is  isotropic,  then  the  cosine  of  the  scattering 
angle  is  given  by  the  expression  CSA  =  RN  for  reflection  into  the 
upper  hemisphere  and  by  the  expression  CSA  =  -RN  for  reflection  into 
the  lower  hemisphere,  where  RN  is  a  random  number. 

When  reflection  is  anisotropic,  the  scattering  angle  CSA  is 
selected  from  the  NRFCOS(NRB)  cosine  values,  RFLCOS(NRB,I) ,  input 
for  equal  intervals  of  the  cumulative  reflection  distributions.  A 
random  number  is  generated  and  the  values  1-1  and  I  determined  such 
that 

(I-l)/NRFCOS(NRB)i  RN  4 I/NRFCOS (NRB) 

The  cosine  of  the  reflection  angle  is  then  calculated  by  inter¬ 
polating  between  the  two  cosine  values  RFLC0S(NRB,I-1)  and  RFLCOS(NRB,I) 
to  give  a  cosine  corresponding  to  the  random  number. 

3.1.8  Printout  of  Scattered  Intensities 

As  discussed  previously  in  this  section  the  scattered  light  in¬ 
tensity  as  calculated  with  the  LITE-I  code  is  printed  in  several 
different  formats.  The  scattered  intensity  is  given  as  a  function 
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of  the  order  of  collision  at  each  receiver  position  so  that  the 
effect  of  multiple  collisions  may  be  analysed.  In  addition,  the 
angular  distribution  of  the  scattered  intensity  is  printed  as  a 
function  of  the  order  of  reflection  from  the  first  reflection  surface 
listed  in  the  problem  input  data.  Thus  the  intensities  listed  for 
the  zero  order  of  reflection  include  the  estimates  of  the  light  in¬ 
tensity  from  all  collisions  that  occur  before  a  light  particle  inter¬ 
sects  the  first  reflection  surface.  The  intensities  listed  for  the 
first  order  of  reflection  include  those  estimates  of  the  light  inten¬ 
sity  reflected  into  a  receiver  due  to  the  first  reflection  and  the 
estimates  of  the  intensity  from  all  those  collisions  that  occur  in 
the  atmosphere  after  the  particle  has  one  reflection  and  before  it 
is  reflected  a  second  time. 

The  printout  of  the  scattered  intensity  as  a  function  of  the 
order  of  reflection  allows  one  to  use  the  results  for  a  problem  .run 
with  one  ground  or  cloud  albedo  to  predict  the  intensity  for  a  prob¬ 
lem  ia  which  all  other  parameters  remain  the  same  except  the  ground 
or  cloud  albedo.  This  allows  a  determination  of  the  effect*  of  vary¬ 
ing  the  ground  or  cloud  albedo  with  only  one  computer  run.  To  con¬ 
vert  the  printout  from  one  albedo  to  another,  the  intensities  listed 
for  the  nth  order  of  reflection  should  be  multiplied  by  the  ratio 
(ALBEDO '/ALBEDO)11 

where  ALBEDO  is  the  albedo  used  to  obtain  the  printout, 

ALBEDO'  is  the  albedo  for  which  the  intensities  are  desired,  and 
n  is  the  order  of  reflection. 
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After  multiplying  the  intensities  listed  for  each  order  of  re¬ 
flection  by  the  above  ratios,  the  products  may  be  summed  over  all 
orders  of  reflection  to  give  the  total  intensity  for  the  problem 
with  the  albedo  ALBEDO'.  If  the  t’-ove  method  is  to  be  used  to  study 
the  effect  of  varying  the  ground  or  cloud  albedo,  it  would  be  advis¬ 
able  to  run  the  initial  problem  for  the  highest  albedo  to  be  con¬ 
sidered. 


The  scattered  intensities  are  also  printed  as  a  function  of  the 
region  of  scatter  for  each  receiver  point  to  determine  those  regions 
of  the  atmosphere  in  which  the  scattering  is  most  significant.  Thus, 
if  a  series  of  similar  problems  are  to  be  run,  the  first  problems 
may  be  used  to  designate  the  more  important  regions  of  scatter  and 
those  regions  may  be  emphasized  in  later  problems  by  assigning  those 
regions  a  higher  importance  number. 

3.1.9  Direct  Intensity  Calculation 


The  direct  intensity  at  each  receiver  point  is  calculated  with 


the  equation 


DI(J)  -  DBSS(J)e 


-(TAUHD(J)  -  TAUHS) /cos6 


where  DI(J)  is  the  direct  intensity  at  the  Jth  receiver  position, 
DBSS(J)  is  the  number  of  light  particles  per  unit  solid  angle 

emitted  from  the  source  in  the  direction  toward  the  Jth 


receiver  per  unit  source  strength, 

TAUHD(J)  is  the  height  of  the  Jth  detector  in  mean-free-path 
lengths, 

TAUHS  is  the  source  height  in  mean-free-path  lengths. 


50 


cos6  is  the  cosine  of  the  angle  between  the  line  joining  the 
source  and  receiver  points  and  the  vertical  vector,  and 
T  is  the  distance  from  the  source  to  the  receiver  point. 

When  the  source  is  isotropic,  DBSS(J)  *  l/4n  for  all  receiver  points. 
If  the  source  is  anisotropic,  then 


DBSS(J) 


PAG(I)  -  PAG(I-l) 

2* |CANG(I)  -  CANG(I-l) | 


where  CANG(I) £  cos© <  CANG(I-l)  and  PAG(I)  is  the  probability  that  a 
particle  will  be  emitted  from  the  source  in  the  polar  angle  interval 
0°  to  cos”1(CANG(I)> . 


3.2  Description  of  the  LITE-II  Code 

The  methods  used  in  the  LITE-II  code  to  treat  the  source  angu¬ 
lar  distribution  and  biasing  of  the  source  angular  distribution  are 
the  same  as  the  methods  previously  described  for  the  LITE-I  code, 
except  that  the  angular  distribution  data  input  to  the  LITE-II  code 
are  for  a  plane  source  rather  than  a  point  source.  Thus  it  is  possi¬ 
ble  to  treat  plane  sources  with  an  arbitrary  angular  distribution 
located  at  the  top  of  the  atmosphere  or  anywhere  withir  the  atmosphere. 

The  methods  used  to  select  random  optical  path  lengths  between 
collisions  in  the  LITE-II  code  are  the  same  as  those  described  in 
Section  3.1  for  the  LITE-I  code.  The  cross  section  input  data  and 
geometry  description  for  the  LITE-II  code  also  correspond  to  that 
used  in  the  LITE-I  code. 


The  major  difference  between  the  two  codes  is  in  the  methods 
used  to  estimate  the  scattered  intensity  at  the  receiver  positions. 
In  the  LITE-II  code,  estimates  of  the  scattered  intensity  at  plane 
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receivers  are  made  rather  than  at  point  receive  Through  the  use 
of  the  reciprocity  principle,  the  estimates  of  t  a  intensity  from  a 
point  source  at  infinite  plane  receivers  are  applied  to  predict  the 
intensity  at  point  receivers  due  to  an  infinite  plane  source. 

In  the  LITE-II  cede  a  new  particle  direction  after  each  colli¬ 
sion  is  chosen  from  either  the  Rayleigh  or  Mie  scattering  phase  func¬ 
tion  as  discussed  in  Section  3.1.2.  Then  an  estimate  of  the  light 
intensity  reaching  each  of  the  receiver  planes  that  Intersect  the 
particle's  positive  directional  vector  is  made.  The  estimate  of  the 
intensity  at  a  receiver  plane  is  given  by  the  equation, 

RESULT  -  (WAIT*e_RH0T) / j  cos 6 1 
where  WAIT  is  the  particle  weight  after  collision, 

RHOT  is  the  number  of  mean-free-path  lengths  between  the  colli¬ 
sion  point  and  receiver  plane  along  the  particle's  direc¬ 
tional  vector,  and 

cos6  is  the  angle  between  the  particle's  directional  vector 
and  the  positive  vertical  vector. 

In  the  LITE-II  code  the  estimates  of  the  intensity  at  each  re¬ 
ceiver  plane  are  stored  and  later  printed  as  a  function  of  the  cosine 
of  the  angle  between  the  particle's  direction  as  it  crosses  the  re¬ 
ceiver  plane  and  the  normal  to  the  plane. 

The  equation  used  to  calculate  the  direct  Intensity  in  the  LITE-II 
code  is  the  same  as  that  used  for  the  LITE-1  code.  However,  the  equa¬ 
tion  can  be  made  to  apply  to  plane  parallel  sources  by  inputting  the 
proper  value  for  DBSS(J).  In  order  for  the  direct-beam  calculation 
given  by  the  equation, 


DI(J)  -  DBSS(J)  e--g-A-^Jj-j..TAUHSL'cQ88  _ 

to  apply  to  a  plane  parallel  source,  the  values  RD(J)  Input  for  the 
radial  positions  of  the  receiver  positions  should  be  determined  by 
the  expression 

RD(J)  -  (HD(J)  -  HS)/cos9 

where  HD(J)  is  the  vertical  position  of  the  receiver  plane, 

HS  is  the  vertical  position  of  the  source,  and 
cos9  is  the  cosine  of  the  angle  that  plane  parallel  source 
makes  with  the  vertical  vector. 


In  addition,  the  values  input  for  DBSS(J)  should  be  the  number 

of  spurce  photons  incident  per  unit  area  upon  the  atmosphere  times 

the  secant  of  the  incident  angle  times  the  square  of  the  slant  dis- 
2 

tance  T  between  the  source  and  receiver  positions,  that  is,  DBSS(J) 
2 

T  sec®  for  a  source  strength  of  one  particle  per  unit  area  on  top  of 
the  atmosphere. 


In  the  special  case  where  the  plane  source  has  an  input  angular 
distribution  defined  by  the  parameters  PAG(J)  and  CANG(J),  then  the 
direct-beam  Intensity  at  the  Ith  detector  plane  per  photon  incident 
per  unit  area  of  the  source  can  not  be  evaluated  by  the  code.  It 
can,  however,  be  computed  separate  from  the  Monte  Carlo  calculations 
by  use  of  the  expression 

DI(I)  -  2"  /  G(e)e"T(I)sec0taned0 
0 


CANG(J) <  cos9  iCANG(J-l)  . 
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The  function  G(6)  defines  the  number  of  source  particles  per  unit 
area  of  the  source  plane  that  are  moving  per  unit  solid  angle  in 
the  direction  0,  where  6  is  the  source  polar  angle  measured  from  the 
vertical  H  axis.  When  the  radiation  incident  on  the  top  of  the  at¬ 
mosphere  is  isotropic,  then  G(6)  •  If  the  plane  source  is  located 
within  the  atmosphere  and  is  isotropic  in  all  directions,  then 

G<e>  -  • 
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IV.  COMPARISONS  WITH  LITE  CODE  RESULTS 

In  order  to  insure  that  the  LITE  code  calculates  properly  the 
intensity  of  scattered  light,  several  test  problems  have  been  run  to 
compare  with  other  calculations.  Calculated,  rather  than  experimen¬ 
tal  data,  were  chosen  to  validate  the  LITE  cc^es,  since  in  most  cases, 
input  information  for  the  calculated  results  are  much  better  defined 
than  are  the  atmospheric  and  terrain  conditions  under  which  the  ex¬ 
perimental  data  were  taken.  Although  some  of  the  calculated  data 
chosen  for  this  purpose  do  not  pertain  to  light  transport  in  the  at¬ 
mosphere,  the  scattering  mechanisms  involved  in  the  calculational 
methods  are  the  same  as  thoce  used  in  the  LITE  codes. 

4.2  LITE-I  Code  Results 

Calculations  were  made  using  the  K74  Monte  Carlo  transport  code 
to  determine  the  direct  and  scattered  intensity  as  a  function  of  dis¬ 
tance  from  a  point  isotropic  source  in  an  infinite  media  for  com¬ 
parison  with  similar  results  obtained  from  the  LITE-I  code. 

Figure  8  shows  the  comparison  of  the  LITE-I  code  results  with 
the  K74  code  results  for  a  .45  micron  wave  length  point  isotropic 
source  embedded  with  an  infinite  air  medium  in  which  the  aerosol 

3 

size  distribution  was  taken  to  be  proportional  to  1/r  ,  and  the  ratio 

of  Rayleigh  to  total  scattering  cross  section  was  taken  to  be  0.58. 

Results  for  problems  run  with  the  ALGOL  version  of  the  LITE-I  code 

are  shown.  The  LITE-I  code  calculated  intensities  for  point  receivers 

2 

located  a  distance  R  from  the  source  were  multiplied  by  4irR  ,  since 
this  was  the  form  of  the  K74  results.  These  data  are  for  ten  orders 
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oz  scattering  which  does  not  represent  the  total  scattered  inten¬ 
sity  in  this  case,  because  in  an  infinite  media,  orders  of  scatter¬ 
ing  above  twenty  would  still  contribute  significantly  to  the  total 
scattered  intensity. 

Several  LITE-I  code  problems  were  run  to  compare  with  Coulson, 
Dave,  and  Sekera's  calculations  of  the  radiation  emerging  from  a 
planetary  atmosphere  with  Rayleigh  scattering  (Reference  10). 

Coulson’ s  data  were  for  plane  monodirectional  sources  incident  upon 
the  top  of  the  atmosphere.  In  order  to  compare  the  LITE-I  code  re¬ 
sults  with  the  plane  source  data,  the  scattered  intensities  I(x,a,p) 
as  a  function  of  the  radial  distance  from  the  vertical  axis  through 
a  point  monodirectional  source  were  calculated,  and  these  intensities 
were  then  integrated  tq  give  the  results  for  a  plane  source.  The 
radial  distribution  of  the  transmitted  scattered  intensity  is  denoted 
as  I(t ,a,p) ,  where  t  is  the  thickness  of  the  atmosphere  in  mean-free- 
paths,  a  is  the  ground  albedo,  and  p  is  the  radial  distance  along  the 
ground  surface  from  the  vertical  axis  through  the  source  point.  The 
integral 

00 

T(x,a)  ■  2ir  J  pl(x,a,p)dp 
0 

gives  the  scattered  intensity  transmitted  through  a  planetary  atmos¬ 
phere  c f  thickness  x  mean-free-paths  and  ground  albedo  of  a  for  a 
plane  source  incident  on  top  of  the  atmosphere.  Also  of  primary 
interest  in  running  the  LITE  problems  to  compare  with  the  data  of 
Coulson  et  al. ,  was  the  effect  of  varying  the  extinction  coefficient 
with  altitude.  To  check  the  effect  of  varying  the  extinction  co¬ 
efficient  with  altitude,  hand  calculations  were  made  of  the  single 
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scattered  intensity  transmitted  through  a  Rayleigh  atmosphere  of 

0*5  optical  thickness  in  which  the  extinction  coefficient  was  held  \ 

constant  with  altitude  and  through  a  Rayleigh  atmosphere  of  the  same 

optical  thickness  in  which  the  extinction  coefficient  was  allowed 

to  vary  exponentially  with  altitude  according  to  the  expression 

E(h)  -  0.00625e°,0125h  , 
where  h  is  in  arbitrary  units. 

The  hand  calculations  are  compared  with  the  single-scattered 
intensities  calculated  with  the  LITE  code  in  Figure  9.  In  both 
instances  the  agreement  is  very  good.  The  variable  density  problem 
overpredicted  the  intensities  for  the  constant  density  problem  for 
detector  points  near  the  vertical  axis  and  underpredicted  the  con¬ 
stant  density  results  for  detector  points  at  the  larger  distances. 

Of  significant  Interest  is  the  fact  that  the  integration  of  the  re¬ 
sults  for  the  two  cases  over  the  lower  plane  of  the  atmosphere  gives 
essentially  the  same  results.  Thus  it  would  indicate  that  in  situations 
involving  the  transmission  of  radiation  through  a  slab  geometry  due 
to  a  plane  source  incident  on  the  slab,  only  the  optical  thickness  of 
the  slab  is  of  significance;  but  in  situations  involving  a  point 
source,  the  manner  in  which  the  extinction  coefficient  varies  with 
altitude  does  affect  the  intensity  at  a  distance,  p,  from  the  verti¬ 
cal  axis  through  the  source  point* 

Two  LITE-I  problems  were  run  for  a  Rayleigh  atmosphere  of  0*13 
optical  thickness.  For  both  problems,  a  monodirectional  point  source 
was  incident  normal  to  an  atmosphere  that  decreased  in  density  with 
increasing  altitude.  In  the  first  problem  the  ground  albedo  was  zero, 
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Fig.  9.  Comparison  of  Hand-Calculated  and  LITE-1  Code  Single  Scattered  Intensities: 
t  =  0.5,  60  =  0°,  Albedo  =  0.0 
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and  in  Che  second  problem  the  ground  albedo  was  0.8.  The  results  of 
the  two  problems  are  shown  in  Figure  10  By  integrating  the  results 
of  the  two  problems  over  a  plane  surface,  the  transmitted  intensities 
for  a  plane  source  incident  normal  to  a  0  13  mean-free-path  thick 
Rayleigh  atmosphere  were  obtained  for  ground  albedos  of  zero  and  0.8. 
These  transmitted  intensities  are  compared  with  the  transmitted  in¬ 
tensities  obtained  from  Coulson's  data  in  Table  II 

Table  II.  Scattered  Light  Transmitted  Through  Rayleigh  Atmospheres 
for  a  Normal  Incident  Plane  Parallel  Broad  Beam  Source 


T 

Air 

Density 

Albedo 

Scattered  Intensity 

-2  -2 
(photon  cm  /source  photon  cm  ) 

LITE  Code 

Coulson 

0.5 

Constant 

0.0 

0.405 

0.403 

0.5 

Variable 

0.0 

0.358 

0.403 

0.5 

Constant 

0.8 

0.964 

1.032 

0.5 

Variable 

0.8 

1.06 

1.032 

0.13 

Variable 

0.0 

0.160 

0.1655 

0.13 

Variable 

0.8 

0.476 

0.440 

Also  shown  in  Table  II  are  the  integrated  intensities  for  two 
problems  run  for  Rayleigh  atmospheres  of  0.5  optical  thickness.  One 
of  the  problems  considered  an  atmosphere  of  constant  density,  and  the 
other,  a  variable  density  atmosphere  in  which  the  extinction  cross 
section  varied  exponentially  with  altitude  according  to  the  expression, 

1(h)  -  0. 00625"° ,0125h  . 

In  both  problems  the  ground  reflection  distribution  was  input 
to  conform  to  Lambert's  law  of  reflection.  The  printout  of  the 
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Fig.  10.  Scattered  Intensity  Transmitted  through  a  0.13  Mean-Free-Path  Thick  Variable  Density  Royleigh 
Atmosphere:  Normal  Incident  Line  Beam  Souice 
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scattered  intensities  as  a  function  of  the  order  of  reflection  for 
a  ground  albedo  of  0.8  made  it  possible  to  determine  the  scattered 
intensities  for  a  ground  albedo  of  0.0c  The  integrated  intensities 
obtained  in  this  manner  for  a  ground  albedo  of  0.0  are  also  shown 
in  Table  II.  The  agreement  between  the  LITE-I  code  data  and  Coulson's 
data  is  reasonably  good  for  both  the  constant  and  variable  density 
cases.  There  is  a  possibility  of  some  error  in  Integrating  the  LITE-I 
code  data,  because  the  spacing  of  the  receiver  points.  The  intensity 
as  a  function  of  radial  position  decreases  very  rapidly  near  the 
source  axis  and  a  slight  shift  of  the  radial  distribution  curve  in 
this  area  can  change  the  integrated  intensity  by  as  much  as  10Z. 

The  sharp  decrease  can  be  observed  in  Figures  11  and  12,  where 
radial  distributions  of  the  transmitted  scattered  Intensities  ob¬ 
tained  with  the  two  problems  are  shown. 

The  difference  in  the  radial  intensity  of  the  multiple  scattered 

radiation  for  the  constant  and  variable  density  atmosphere  is  seen 

to  follow  the  same  pattern  as  did  the  radial  intensities  for  the 

single  scattered  intensities.  In  each  of  the  two  problems  plotted 

in  Figures  11  and  12,  a  maximum  of  20  collisions  per  history  was 

allowed.  The  contributions  from  the  20th  order  of  collision  were 

3  5 

a  factor  of  approximately  10  to  10  lower  than  those  from  the  first 
order  of  collision.  Slightly  more  than  one-half  of  the  particles 
escaped  out  of  the  problem  geometry  before  they  had  undergone  20 
collisions. 

Problems  have  also  been  run  with  the  LITE-I  code  to  produce 
results  for  comparison  with  exact  solutions  reported  by  Beach 
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RADIAL  DISTANCE  FROM  VERTICAL  AXIS  (arbitrary  units) 


Fig.  II.  Scattered  Intensity  Transmitted  through  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmosphere:  Normal 
Incident  Line  Beam  Source,  Albedo  0.8 
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RADIAL  DISTANCE  FROM  VERTICAL  AXIS  (arbitrary  units) 


Fig.  12.  Scattered  Intensity  Transmitted  through  a  0.5  Mean-Free-Poth  Thick  Rayleigh  Atmosphere:  Normal 
Incident  Line  Beam  Source,  Albedo  -  0.0 
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et  al.  (Reference  18)  for  one-velocity  neutron  diffusion  problems. 
The  first  problem  considered  was  for  a  plane  isotropic  source  em¬ 
bedded  within  an  infinite  medium  in  which  scattering  was  isotropic 
in  the  laboratory  system,  and  the  scattering- to- total  cross  sec¬ 
tion  ratio  was  0.5.  Beach  reported  calculations  of  the  total  in¬ 
tensity  as  a  function  of  distance  from  the  plane  source.  The  LITE-I 
code  was  run  to  give  the  intensity  at  point  receivers  located  at 

several  distances  from  a  point  isotropic  source.  The  intensity  was 

2 

then  multiplied  by  4xR  and  plotted  as  a  function  of  the  distance 
R  from  the  point  source  as  shown  in  Figure  13.  It  is  thought  that 
the  point  at  900  units  is  probably  an  underprediction  of  the  inten¬ 
sity,  since  sampling  at  this  distance  from  the  source  would  be  in¬ 
adequate  unless  some  biasing  scheme  is  used.  Therefore,  the  curve, 
as  drawn  through  the  previous  four  points,  was  extrapolated  to  pre¬ 
dict  the  intensity  at  900  units  and  larger  distances. 

To  apply  the  LITE-I  code  results  to  obtain  the  intensity  at  a 
distance  h  below  a  plane  isotropic  source,  the  following  integral 


I(h)  -  /  2ffI0G(R)PdP 

0  4ttR  2 


4.1.1 


An  DISTANCE  SQUARED  TIMES  INTENSITY  (4irparticles/source  particle) 
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Fig.  13.  4rrR2  Times  Intensity  as  o  Function  of  Distance  from  o  Point  Isotropic  Source  in  on  Infinite  Medium? 
Scattering-to-Tatol  Cross  Section  Ratio  =  0.5 
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where  G(R)  is  4irR  times  the  intensity  at  distance  R  from  a  point 
isotropic  source .  By  making  the  substitutions 

p  ■  h  tan9 
2 

dp  ■  h  sec  6d6 
R  *  h  sec© 

equation  4.1.1  becomes 
T  */2 

1(h)  -  ±2  /  G(R  -  hsec0)tan0d0  4.1.2 

2  0 

Equation  4.1.2  was  Integrated  numerically  for  values  of  h 
corresponding  to  optical  distances  of  t  ■  1.1,  2.2,  4.1,  5.5  and 
8.3.  A  comparison  between  Beach's  data  and  the  integrated  LITE-I 
code  results  is  shown  in  Figure  14. 

A  second  comparison  was  made  between  LITE  code  results  and 
Beach '8  data  for  a  plane  parallel  source  embedded  within  an  in¬ 
finite  media  in  which  all  scattering  was  taken  to  be  isotropic, 
and  the  scattering-to-total  cross  section  ratio  was  0.5.  The  LITE-I 
code  was  used  to  calculate  the  scattered  intensity  at  points  on 
three  planes  located  1,  2  and  3  optical  distances  below  a  point 
monodirectional  source  directed  normal  to  the  receiver  planes.  An 
integration  of  the  intensities  over  the  points  in  each  plane  was 
performed  to  obtain  the  total  scattered  intensity  received  by  each 
plane.  These  integrated  intensities  were  related  to  the  intensity 
from  a  plane  monodirectional  source  through  the  use  of  the  recipro¬ 
city  principle.  A  comparison  of  the  scattered  intensities  derived 
from  the  LITE  code  results  and  those  computed  by  Beach  is  shown  in 
Figure  15. 


INTENSITY  (photons  cm“2/2  photons  cm” 
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Fig.  14.  Intensity  Versus  Distance  from  a  Plane  Isotropic  Source  in  an  Infinite  Medium:  Scattering-ta-Tatal 
Cross  Section  Ratio  =  0.5 


DISTANCE  FROM  SOURCE  PLANE  (meon-free-poths) 


Fig.  15.  Scattered  Intensity  as  a  Function  of  Distonce  from  a  Pione  Porollel  Monodirectional  Source  in  on 
Infinite  Medium:  0Q  =  0° ,  Scattering-to-Totol  Cross  Section  Ratio  =  0.5 
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The  good  comparison  between  Beach's  data  and  the  LITE-I  code 
indicates  that  the  routines  for  treating  Mie  scattering  are  opera¬ 
ting  properly ,  since  the  isotropic  scattering  distributions  vrere 
described  in  the  problem  input  as  Mie  scattering. 

4.?  LITE-I I  Code  Results 

Two  problems  for  a  0.5  mean-free-path  thick  Rayleigh  atmosphere 
were  ruu  with  the  LITE-II  code  to  obtain  results  to  compare  with 
Coulson's  data.  The  first  problem  vras  the  same  as  the  variable 
density  problem  run  with  the  LITE-I  code  for  a  normal  incident 
plane  source  and  a  0.8  ground  albedo.  The  second  problem  was  for  a 
plane  source  incident  upon  a  0.5  mean-free-path  thick  variable  den¬ 
sity  Raleigh  atmosphere  at  an  angle  of  0q  *  cos  ^0.^.  The  second 
problem  was  also  run  for  a  ground  albedo  of  0.8.  The  results  of 
the  LITE-II  code  gives  directly  the  intensity  for  a  plane  source 
incident  upon  the  atmosphere,  so  no  integration  of  the  output  data 
was  necessary.  For  the  two  LITE-II  problems,  one  receiver  was 
placed  in  the  atmosphere  just  slightly  above  the  ground  surface  to 
record  the  transmitted  scattered  Intensity.  A  second  receiver  was 
located  just  slightly  below  the  surface  defining  the  upper  limit 
of  the  atmosphere  to  record  the  Intensity  reflected  back  out  of  the 
top  of  the  atmosphere.  Although  the  two  problems  were  run  for  a 
ground  albedo  of  0.8,  the  results  were  also  converted  to  data  for 
ground  albedos  of  0.0  and  0.25  tor  comparison  with  Coulson's  data 
for  those  albedos.  Tables  III  and  IV  show  a  comparison  of  the  data 
from  the  LITE-II  code  problems  with  Coulson's  data.  The  data  tabu¬ 
lated  in  Tables  III  and  IV  give  the  scattered  intensity  transmitted 
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Table  III.  Scattered  Radiation  Emerging  from  Upper  and  Lower 

Surfaces  of  a  0.5  Mean-Free  Path  Thick  Rayleigh  At¬ 
mosphere  Due  to  a  Normal  Incident  Plane  Source 


Upper  Surface 

Lower 

Surface 

Albedo 

LITE  Code 

Coulson  Data 

LITE  Code 

Coulson  Data 

0.0 

.4442 

.4422 

.3694 

.4025 

0.25 

.7007 

.7106 

.532 

.5672 

0.80 

1.458 

1.486 

.9944 

1.0321 

Table  IV. 

Scattered  Radiation  Emerging  from  the  Upper  and  Lower 

Surfaces  of  a  0.5  Mean-Free-Path  Thick  Rayleigh  Atmo¬ 
sphere  Due  to  Plane  Source  Incident  at  0q  *  cos”^0.6. 

A  1  — 

Upper  Surface 

Lower 

Surface 

Albedo 

LITE  Code  Coulson  Data 

LITE  Code 

Coulson  Data 

0.0 

.43  .4344 

.38 

.3768 

0.25 

.57  .5900 

.47 

.4424 

0.80 

.98  1.0028 

.73 

.7132 
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through  or  reflected  from  the  atmosphere  due  to  a  unit  intensity 
incident  ion  the  top  of  the  atmosphere.  Angular  distributions  of 
the  LITE-II  calculated  scattered  intensities  emerging  from  the  at¬ 
mosphere  were  also  compared  with  Coulson's  data.  Figures  16  through 
27  present  comparisons  of  the  LITE-11  code  calculated  angular  dis¬ 
tributions  at  the  upper  and  lower  surfaces  of  the  atmosphere  for 
each  of  the  three  ground  albedos  for  the  normal  incident  plane  source 
and  the  source  incident  at  6q  •  cos  ^0.6.  Coulson's  angular  distri¬ 
butions  are  given  for  a  source  term  of  ir  photons  per  unit  area  normal 
to  the  direction  of  travel.  The  angular  distribution  data  obtained 
from  the  LITE-II  output,  have  been  multiplied  by  the  same  source 
term  to  produce  the  comparisons  shown  in  Figures  16  through  27.  The 
curves  represent  the  intensity  per  steradian  emitted  from  the  at¬ 
mosphere  as  a  function  of  the  cosine  of  the  angle  measured  from  the 
surface  normal. 

Two  C18  problems  were  run  to  provide  data  on  the  light  inten¬ 
sity  as  a  function  of  depth  in  a  Rayleigh  atmosphere  for  comparison 
with  the  results  of  a  LITE-II  code  problem.  The  thickness  of  the 
atmosphere  was  taken  to  be  0.5  mean-free-paths.  The  ground  albedo 
was  taken  to  be  zero  and  the  source  was  taken  to  be  inc-1 *ent  normal 
to  the  top  of  the  atmosphere.  In  the  first  C18  problem  the  atmo¬ 
spheric  density  was  taken  to  be  constant  with  altitude.  In  the 
second  C18  problem,  the  slab  atmosphere  was  divided  into  5  subslabs 
whose  thicknesses  were  defined  by  the  equation 

0.1  -  /*  0.00625  e°-0125hdh,  1  -  1,5 
hi-l 

h^_i  ■  0  for  i  •  1, 
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Fig.  17.  Angular  DUtribution  oi  the  Scattered  Radiation  Tranimitted  through  a  0.5  Mean-Free-Peth  Thick  Rayleigh  At  men  ph  ere: 
%=&>,  Ground  Albedo  =  0.25 


o 


Fig.  18.  Angulor  Distribution  of  the  Scattered  Rodiation  Transmitted  through  o  0.5  Meon- Free- Path  Thick  Royleigh  Atmosphere: 
0°  ,  Ground  Albedo  “0.8 
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Fig.  19.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5  Meon- Free -Path  Thick  Rayleigh  Atmosphere: 
0o=O°  Ground  Albedo  =  0.0 


o 


Fig.  20.  Angular  Distribution  of  the  Rodiation  Reflected  from  a  0.5  Mean -Free -Path  Thick  Rayleigh  Atmosphere: 
0o~O° ,  Ground  Albedo  -  0.25 


77 


Fig.  21.  Angular  Diitribution  af  the  Radiation  Reflected  from  a  0.5  Mean -Free -Path  Thick  Rayleigh  Atmos| 
6^-0°,  Ground  Albeda  =  0.8 
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Angular  Distribution  of  the  Scattered  Radiation  Transmitted  through  a  0.5  Mean -Free -Path  Thick  Rayleigh  Atmosphere: 
0O=  Cos"1 0.6,  Ground  Albedo  =  0.0 
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Fig.  23.  Angular  Distribution  of  the  Scattered  Radiation  Transmitted  through  a  0.5  Mean -Free -Path  Thick  Rayleigh  Atmosphere: 
Cos“l0.6#  Ground  Albedo  ■  0.25 
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Fig.  24.  Angular  Distribution  of  the  Scattered  Rodiation  Transmitted  through  a  0.5  Mean -Free- Path  Thick  Rayleigh  Atmosphere: 
0Q  =  Cos~’  0.6,  Ground  Albedo  =0.8 
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Fig.  25.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5  Mean -Free -Path  Thick  Rayleigh  Atmosphere: 
Oq-  Cos'1 0.6,  Ground  Albedo  =  0.0 
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Fig.  27.  Angular  Distribution  of  the  Radiation  Reflected  from  a  0.5  Mean- Free -Pa  *h  Thick  Rayleigh  Atmosphere: 
6q~  Cos-' 0.6,  Ground  Albedo  =0.8 
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where  h^_^  and  are  the  distances  from  the  bottom  of  the  atmosphere 
to  the  bottom  and  top  of  each  of  the  five  subslabs  respectively. 

The  results  from  the  two  C18  problems  are  presented  in  Figure  28 
as  a  function  of  the  detector  depth  into  the  atmosphere  in  mean-free- 
paths.  The  results  from  a  l.ITE-II  code  problem  in  which  the  atmo¬ 
spheric  density  was  assumed  to  vary  exponentially  with  altitude  are 
also  presented  in  Figure  28.  It  is  seen  that,  within  the  statisti¬ 
cal  variation  of  the  Monte  Carlo  results,  there  is  no  significant 
differences  between  the  calculated  spatial  distributions  of  the  total 
intensity  as  given  by  the  two  different  codes.  When  the  data  in 
Figure  28  were  plotted  as  a  function  of  distance  in  units  of  h  in¬ 
stead  of  in  mean-free-paths,  it  was  found  that  the  intensity  at  a 
given  distance  is  dependent  on  the  atmospheric  density  distribution 
assumed  for  each  problem.  The  results  of  these  comparisons  also 
show  that  the  emergent  radiation  intensity  at  the  top  and  bottom  of 
the  atmosphere  is  independent  of  the  atmospheric  density  distribu¬ 
tion  assumed  for  each  problem. 

The  angular  distribution  of  the  scattered  intensity  emerging 
from  the  bottom  of  a  0.5  mean-free-path  thick  Rayleigh  atmosphere 
as  given  by  the  C18  calculations  for  a  variable  air  density,  are 
presented  in  Figure  16,  where  they  can  be  compared  with  similar  re¬ 
sults  obtained  from  the  LITE-II  code  and  from  Coulson’s  data.  It 
is  seen  that  results  from  the  three  different  calculations  are  in 
reasonably  good  agreement  considering  the  small  sample  sizes  used 
in  the  LITE-II  and  C18  Monte  Carlo  calculations. 


INTENSITY  (photons  cm~Vphoton  cm 
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DEPTH  INTO  ATMOSPHERE  (mean-free-paths) 


Fig.  28.  Spatial  Distribution  af  the  Radiation  Intensity  as  a  Function  af  Depth  into  a  0.5  Mean -Free -Path  Thick 
Rayleigh  Atmosphere*  Normal  Incident  Plane  Parallel  Source 
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